diff --git a/src/utils_trust_region/algo_trust.irp.f b/src/utils_trust_region/algo_trust.irp.f
index bfb3a0ff..eac17275 100644
--- a/src/utils_trust_region/algo_trust.irp.f
+++ b/src/utils_trust_region/algo_trust.irp.f
@@ -137,9 +137,9 @@ subroutine trust_region_step_w_expected_e(n,H,W,e_val,v_grad,prev_criterion,rho,
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the step and the expected criterion/energy after the step
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -210,9 +210,9 @@ subroutine trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, cri
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute if the step should be cancelled
-  !END_DOC
+  END_DOC
 
   implicit none
  
diff --git a/src/utils_trust_region/apply_mo_rotation.irp.f b/src/utils_trust_region/apply_mo_rotation.irp.f
index a313769d..42aac5f7 100644
--- a/src/utils_trust_region/apply_mo_rotation.irp.f
+++ b/src/utils_trust_region/apply_mo_rotation.irp.f
@@ -25,9 +25,9 @@ subroutine apply_mo_rotation(R,prev_mos)
   
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the new MOs knowing the rotation matrix
-  !END_DOC
+  END_DOC
 
   implicit none
 
diff --git a/src/utils_trust_region/rotation_matrix.irp.f b/src/utils_trust_region/rotation_matrix.irp.f
index 0cf47355..4738fd67 100644
--- a/src/utils_trust_region/rotation_matrix.irp.f
+++ b/src/utils_trust_region/rotation_matrix.irp.f
@@ -61,10 +61,10 @@ subroutine rotation_matrix(A,LDA,R,LDR,n,info,enforce_step_cancellation)
 
   implicit none
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Rotation matrix to rotate the molecular orbitals.
   ! If the rotation is too large the transformation is not unitary and must be cancelled.
-  !END_DOC
+  END_DOC
 
   include 'pi.h'
 
diff --git a/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f b/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f
index 75d04352..bdd1f6ba 100644
--- a/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f
+++ b/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f
@@ -32,9 +32,9 @@
 
 subroutine sub_to_full_rotation_matrix(m,tmp_list,tmp_R,R)
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the full rotation matrix from a smaller one
-  !END_DOC
+  END_DOC
   
   implicit none
 
diff --git a/src/utils_trust_region/trust_region_expected_e.irp.f b/src/utils_trust_region/trust_region_expected_e.irp.f
index 1a18dfb3..b7d849d1 100644
--- a/src/utils_trust_region/trust_region_expected_e.irp.f
+++ b/src/utils_trust_region/trust_region_expected_e.irp.f
@@ -33,9 +33,9 @@ subroutine trust_region_expected_e(n,v_grad,H,x,prev_energy,e_model)
    
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the expected criterion/energy after the application of the step x
-  !END_DOC
+  END_DOC
 
   implicit none
 
diff --git a/src/utils_trust_region/trust_region_optimal_lambda.irp.f b/src/utils_trust_region/trust_region_optimal_lambda.irp.f
index 71f6381d..f007a860 100644
--- a/src/utils_trust_region/trust_region_optimal_lambda.irp.f
+++ b/src/utils_trust_region/trust_region_optimal_lambda.irp.f
@@ -153,9 +153,9 @@ subroutine trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Research the optimal lambda to constrain the step size in the trust region
-  !END_DOC
+  END_DOC
 
   implicit none
   
@@ -508,9 +508,9 @@ function d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta)
   use omp_lib
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -644,9 +644,9 @@ function d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta)
   use omp_lib
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2
-  !END_DOC
+  END_DOC
   
   implicit none
 
@@ -792,9 +792,9 @@ function f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute ||x(lambda)||^2
-  !END_DOC
+  END_DOC
   
   implicit none
 
@@ -903,9 +903,9 @@ function d1_norm_trust_region(n,e_val,w,v_grad,lambda,delta)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 
-  !END_DOC
+  END_DOC
   
   implicit none
 
@@ -1000,9 +1000,9 @@ function d2_norm_trust_region(n,e_val,w,v_grad,lambda,delta)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -1102,9 +1102,9 @@ function f_norm_trust_region(n,e_val,tmp_wtg,lambda)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute ||x(lambda)||^2
-  !END_DOC
+  END_DOC
   
   implicit none
 
@@ -1184,9 +1184,9 @@ function d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta)
   use omp_lib
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -1340,9 +1340,9 @@ function d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta)
   use omp_lib
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -1498,9 +1498,9 @@ function d1_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2
-  !END_DOC
+  END_DOC
 
   implicit none
 
@@ -1592,9 +1592,9 @@ function d2_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2
-  !END_DOC
+  END_DOC
 
   implicit none
 
diff --git a/src/utils_trust_region/trust_region_rho.irp.f b/src/utils_trust_region/trust_region_rho.irp.f
index eb280ef0..45738736 100644
--- a/src/utils_trust_region/trust_region_rho.irp.f
+++ b/src/utils_trust_region/trust_region_rho.irp.f
@@ -47,9 +47,9 @@ subroutine trust_region_rho(prev_energy, energy,e_model,rho)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compute rho, the agreement between the predicted criterion/energy and the real one
-  !END_DOC
+  END_DOC
 
   implicit none
    
diff --git a/src/utils_trust_region/trust_region_step.irp.f b/src/utils_trust_region/trust_region_step.irp.f
index afa19dec..42aa6ed4 100644
--- a/src/utils_trust_region/trust_region_step.irp.f
+++ b/src/utils_trust_region/trust_region_step.irp.f
@@ -375,9 +375,9 @@ subroutine trust_region_step(n,nb_iter,v_grad,rho,e_val,w,x,delta)
 
   include 'pi.h'
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Compuet the step in the trust region
-  !END_DOC
+  END_DOC
 
   implicit none
 
diff --git a/src/utils_trust_region/vec_to_mat_v2.irp.f b/src/utils_trust_region/vec_to_mat_v2.irp.f
index e184d3ba..9140b8d3 100644
--- a/src/utils_trust_region/vec_to_mat_v2.irp.f
+++ b/src/utils_trust_region/vec_to_mat_v2.irp.f
@@ -8,9 +8,9 @@
 
 subroutine vec_to_mat_v2(n,m,v_x,m_x)
 
-  !BEGIN_DOC
+  BEGIN_DOC
   ! Vector to antisymmetric matrix
-  !END_DOC
+  END_DOC
   
   implicit none