diff --git a/src/casscf_cipsi/EZFIO.cfg b/src/casscf_cipsi/EZFIO.cfg
index 18e0b6b1..5b72d906 100644
--- a/src/casscf_cipsi/EZFIO.cfg
+++ b/src/casscf_cipsi/EZFIO.cfg
@@ -79,3 +79,9 @@ type: logical
 doc: If |true|, the pt2_max value in the CIPSI is set to 10-10 and will not change
 interface: ezfio,provider,ocaml
 default: False
+
+[act_mos_opt]
+type: logical
+doc: If |true|, the active orbitals are also optimized variationally
+interface: ezfio,provider,ocaml
+default: False
diff --git a/src/casscf_cipsi/NEED b/src/casscf_cipsi/NEED
index 11d1a78c..32f5ae90 100644
--- a/src/casscf_cipsi/NEED
+++ b/src/casscf_cipsi/NEED
@@ -3,4 +3,4 @@ selectors_full
 generators_cas
 two_body_rdm
 dav_general_mat
-mo_optimization
+mo_optimization_utils
diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f
index dc3e2245..b64a9d8f 100644
--- a/src/casscf_cipsi/casscf.irp.f
+++ b/src/casscf_cipsi/casscf.irp.f
@@ -46,94 +46,101 @@ subroutine run
   do while (.not.converged)
     print*,'pt2_max = ',pt2_max
     call run_stochastic_cipsi(Ev,PT2)
-    print*,'Ev,PT2',Ev(1),PT2(1)
-    E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
-    energy_old = energy
-    energy = eone+etwo+ecore
-    pt2_max_before = pt2_max
-
-    call write_time(6)
-    call write_int(6,iteration,'CAS-SCF iteration = ')
-    call write_double(6,energy,'State-average CAS-SCF energy = ')
-!    if(n_states == 1)then
-!     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
-!     call ezfio_get_casscf_cipsi_energy(PT2)
-     double precision :: delta_E_istate, e_av
-     e_av = 0.d0
-     do istate=1,N_states
-      e_av += state_average_weight(istate) * Ev(istate) 
-      if(istate.gt.1)then
-       delta_E_istate = E_PT2(istate)  - E_PT2(1)
-      write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2    = ',delta_E_istate
-      endif
-      write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate)
-      write(*,'(A6,I2,A18,F16.10)')'state ',istate,'     PT2 energy = ',PT2(istate)
-!      call write_double(6,E_PT2(istate),'E + PT2 energy = ')
-!      call write_double(6,PT2(istate),'  PT2          = ')
-     enddo
-     call write_double(6,e_av,'State-average CAS-SCF energy bis = ')
-     call write_double(6,pt2_max,' PT2_MAX       = ')
+!    if(act_mos_opt)then DOES NOT WORK
+!     call run_orb_opt_trust_v2
+!     call run_stochastic_cipsi(Ev,PT2)
 !    endif
-
-    print*,''
-    call write_double(6,norm_grad_vec2,'Norm of gradients = ')
-    call write_double(6,norm_grad_vec2_tab(1), '     Core-active  gradients = ')
-    call write_double(6,norm_grad_vec2_tab(2), '     Core-virtual gradients = ')
-    call write_double(6,norm_grad_vec2_tab(3), '   Active-virtual gradients = ')
-    print*,''
-    call write_double(6,energy_improvement, 'Predicted energy improvement = ')
-
-    if(criterion_casscf == "energy")then
-     converged = dabs(energy_improvement) < thresh_scf
-    else if (criterion_casscf == "gradients")then
-     converged = norm_grad_vec2 < thresh_scf
-    else if (criterion_casscf == "e_pt2")then
-     delta_E = 0.d0
-     do istate = 1, N_states
-      delta_E += dabs(E_PT2(istate) - ept2_before(istate))
-     enddo
-     converged = dabs(delta_E) < thresh_casscf
-    endif
-    ept2_before = E_PT2
-    if(.not.small_active_space)then
-     if(adaptive_pt2_max)then
-      pt2_max = dabs(energy_improvement / (pt2_relative_error))
-      pt2_max = min(pt2_max, pt2_max_before)
-      if(n_act_orb.ge.n_big_act_orb)then
-       pt2_max = max(pt2_max,pt2_min_casscf)
-      endif
+    if(.True.)then
+     print*,'Ev,PT2',Ev(1),PT2(1)
+     E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
+     energy_old = energy
+     energy = eone+etwo+ecore
+     pt2_max_before = pt2_max
+ 
+     call write_time(6)
+     call write_int(6,iteration,'CAS-SCF iteration = ')
+     call write_double(6,energy,'State-average CAS-SCF energy = ')
+!!    if(n_states == 1)then
+!!     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
+!!     call ezfio_get_casscf_cipsi_energy(PT2)
+      double precision :: delta_E_istate, e_av
+      e_av = 0.d0
+      do istate=1,N_states
+       e_av += state_average_weight(istate) * Ev(istate) 
+       if(istate.gt.1)then
+        delta_E_istate = E_PT2(istate)  - E_PT2(1)
+       write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2    = ',delta_E_istate
+       endif
+       write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate)
+       write(*,'(A6,I2,A18,F16.10)')'state ',istate,'     PT2 energy = ',PT2(istate)
+!!      call write_double(6,E_PT2(istate),'E + PT2 energy = ')
+!!      call write_double(6,PT2(istate),'  PT2          = ')
+      enddo
+      call write_double(6,e_av,'State-average CAS-SCF energy bis = ')
+      call write_double(6,pt2_max,' PT2_MAX       = ')
+!!    endif
+ 
+     print*,''
+     call write_double(6,norm_grad_vec2,'Norm of gradients = ')
+     call write_double(6,norm_grad_vec2_tab(1), '     Core-active  gradients = ')
+     call write_double(6,norm_grad_vec2_tab(2), '     Core-virtual gradients = ')
+     call write_double(6,norm_grad_vec2_tab(3), '   Active-virtual gradients = ')
+     print*,''
+     call write_double(6,energy_improvement, 'Predicted energy improvement = ')
+ 
+     if(criterion_casscf == "energy")then
+      converged = dabs(energy_improvement) < thresh_scf
+     else if (criterion_casscf == "gradients")then
+      converged = norm_grad_vec2 < thresh_scf
+     else if (criterion_casscf == "e_pt2")then
+      delta_E = 0.d0
+      do istate = 1, N_states
+       delta_E += dabs(E_PT2(istate) - ept2_before(istate))
+      enddo
+      converged = dabs(delta_E) < thresh_casscf
      endif
-    endif
-    print*,''
-    call write_double(6,pt2_max, 'PT2_MAX for next iteration = ')
-
-    mo_coef = NewOrbs
-    mo_occ  = occnum
-    if(.not.converged)then
-     call save_mos
-     iteration += 1
-     if(norm_grad_vec2.gt.0.01d0)then
-      N_det = N_states
-     else
-      N_det = max(N_det/8 ,N_states)
-     endif
-     psi_det = psi_det_sorted
-     psi_coef = psi_coef_sorted
-     read_wf = .True.
-     call clear_mo_map
-     SOFT_TOUCH mo_coef N_det psi_det psi_coef
+     ept2_before = E_PT2
      if(.not.small_active_space)then
       if(adaptive_pt2_max)then
-        SOFT_TOUCH pt2_max  
+       pt2_max = dabs(energy_improvement / (pt2_relative_error))
+       pt2_max = min(pt2_max, pt2_max_before)
+       if(n_act_orb.ge.n_big_act_orb)then
+        pt2_max = max(pt2_max,pt2_min_casscf)
+       endif
       endif
      endif
-     if(iteration .gt. 3)then
-      state_following_casscf = state_following_casscf_cipsi_save
-      soft_touch state_following_casscf
+     print*,''
+     call write_double(6,pt2_max, 'PT2_MAX for next iteration = ')
+ 
+     mo_coef = NewOrbs
+     mo_occ  = occnum
+     if(.not.converged)then
+      call save_mos
+      iteration += 1
+      if(norm_grad_vec2.gt.0.01d0)then
+       N_det = N_states
+      else
+       N_det = max(N_det/8 ,N_states)
+      endif
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      read_wf = .True.
+      call clear_mo_map
+      SOFT_TOUCH mo_coef N_det psi_det psi_coef
+      if(.not.small_active_space)then
+       if(adaptive_pt2_max)then
+         SOFT_TOUCH pt2_max  
+       endif
+      endif
+      if(iteration .gt. 3)then
+       state_following_casscf = state_following_casscf_cipsi_save
+       soft_touch state_following_casscf
+      endif
      endif
     endif
-
+ 
   enddo
+  if(.True.)then
      integer :: i
     print*,'Converged CASSCF '
     print*,'--------------------------'
@@ -153,6 +160,7 @@ subroutine run
 !   write(*,*)mcscf_fock_alpha_mo(i,i)
   enddo
 
+ endif
 
 end
 
diff --git a/src/mo_optimization_utils/class.irp.f b/src/mo_optimization/class.irp.f
similarity index 100%
rename from src/mo_optimization_utils/class.irp.f
rename to src/mo_optimization/class.irp.f