From 248545e558effff89eb124b999baab562bb8db8d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 31 May 2021 14:01:33 +0200
Subject: [PATCH] DIR -> OMP

---
 configure.sh  |   4 +-
 src/det.irp.f | 340 +++++++++++++++++++++++++-------------------------
 src/mo.irp.f  | 170 ++++++++++++-------------
 3 files changed, 257 insertions(+), 257 deletions(-)

diff --git a/configure.sh b/configure.sh
index 535eea0..c48dfe7 100755
--- a/configure.sh
+++ b/configure.sh
@@ -39,11 +39,11 @@ export PATH="\${QMCCHEM_PATH}/bin:\${PATH}"
 export LD_LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LD_LIBRARY_PATH}"
 export LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LIBRARY_PATH}"
 export QMCCHEM_MPIRUN="mpirun"
-export QMCCHEM_MPIRUN_FLAGS="--bind-to-core"
+export QMCCHEM_MPIRUN_FLAGS=""
 export C_INCLUDE_PATH="\${QMCCHEM_PATH}/include:\${C_INCLUDE_PATH}"
 #export QMCCHEM_NIC=ib0
 source \${QMCCHEM_PATH}/irpf90/bin/irpman
 #source \${QMCCHEM_PATH}/EZFIO/Bash/ezfio.sh
-eval $(opam env)
+eval \$(opam env)
 EOF
 
diff --git a/src/det.irp.f b/src/det.irp.f
index a7c0fca..91e2034 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1,32 +1,32 @@
 BEGIN_PROVIDER  [ integer, det_i ]
-  
+
   BEGIN_DOC
   ! Current running alpha determinant
   END_DOC
   det_i=det_alpha_order(1)
-  
+
 END_PROVIDER
 
 BEGIN_PROVIDER  [ integer, det_j ]
-  
+
   BEGIN_DOC
   ! Current running beta determinant
   END_DOC
   det_j=det_beta_order(1)
-  
+
 END_PROVIDER
 
 subroutine det_update(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m(LDS) ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,n)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   if (d == 0.d0) then
     return
   endif
@@ -193,15 +193,15 @@ end
 
 subroutine det_update2(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m(2)   ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,2)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,2) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   S(1,l) = m(1)
   S(2,l) = m(2)
   S_inv(1,1) = S(1,1)
@@ -209,37 +209,37 @@ subroutine det_update2(n,LDS,m,l,S,S_inv,d)
   S_inv(2,1) = S(1,2)
   S_inv(2,2) = S(2,2)
   call invert2(S_inv,LDS,n,d)
-  
+
 end
 
 subroutine det_update1(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m(1)   ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,1)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,1) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   S(1,l) = m(1)
   S_inv(1,1) = S(1,1)
   call invert1(S_inv,LDS,n,d)
-  
+
 end
 
 subroutine det_update3(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m(3)   ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,3)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,3) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   integer                        :: i
   do i=1,3
     S(i,l) = m(i)
@@ -249,22 +249,22 @@ subroutine det_update3(n,LDS,m,l,S,S_inv,d)
     S_inv(2,i) = S(i,2)
     S_inv(3,i) = S(i,3)
   enddo
-  
+
   call invert3(S_inv,LDS,n,d)
-  
+
 end
 
 subroutine det_update4(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m(4)   ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,4)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,4) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   double precision               :: u(4), z(4), w(4), lambda, d_inv
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
   integer                        :: i,j
@@ -276,7 +276,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d)
   z(2) = S_inv(1,2)*u(1) + S_inv(2,2)*u(2) + S_inv(3,2)*u(3) + S_inv(4,2)*u(4)
   z(3) = S_inv(1,3)*u(1) + S_inv(2,3)*u(2) + S_inv(3,3)*u(3) + S_inv(4,3)*u(4)
   z(4) = S_inv(1,4)*u(1) + S_inv(2,4)*u(2) + S_inv(3,4)*u(3) + S_inv(4,4)*u(4)
-  
+
   d_inv = 1.d0/d
   d = d+z(l)
   lambda = d_inv*d
@@ -284,7 +284,7 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d)
     d = 0.d0
     return
   endif
-  
+
   !DIR$ VECTOR ALIGNED
   do i=1,4
     w(i) = S_inv(i,l)*d_inv
@@ -297,22 +297,22 @@ subroutine det_update4(n,LDS,m,l,S,S_inv,d)
     S_inv(j,i) = S_inv(j,i)*lambda -z(i)*w(j)
    enddo
   enddo
-  
+
 end
 
 BEGIN_TEMPLATE
 ! Version for mod(n,4) = 0
 subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m($n)  ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,$n)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   double precision               :: u($n), z($n), w($n), lambda, d_inv
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
   !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN
@@ -321,20 +321,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (LDS >= $n)
   integer                        :: i,j
   double precision :: zj, zj1, zj2, zj3
-  
+
   !DIR$ NOPREFETCH
-  !DIR$ SIMD NOVECREMAINDER
+  !OMP$ SIMD NOVECREMAINDER
   do i=1,$n
     u(i) = m(i) - S(i,l)
   enddo
-  
+
   zj = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER
+  !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER
   do i=1,$n-1,4
     zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1)  &
-            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) 
+            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3)
   enddo
 
   d_inv = 1.d0/d
@@ -344,7 +344,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
     d = 0.d0
     return
   endif
-  
+
   !DIR$ VECTOR ALIGNED
   do j=1,$n,4
    zj  = 0.d0
@@ -353,26 +353,26 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj3 = 0.d0
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
+   !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
    do i=1,$n
     zj  = zj  + S_inv(i,j  )*u(i)
     zj1 = zj1 + S_inv(i,j+1)*u(i)
     zj2 = zj2 + S_inv(i,j+2)*u(i)
     zj3 = zj3 + S_inv(i,j+3)*u(i)
    enddo
-   z(j  ) = zj  
-   z(j+1) = zj1 
-   z(j+2) = zj2 
+   z(j  ) = zj
+   z(j+1) = zj1
+   z(j+2) = zj2
    z(j+3) = zj3
   enddo
 
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
+  !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
   do i=1,$n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
   enddo
-  
+
   do i=1,$n,4
    zj  = z(i  )
    zj1 = z(i+1)
@@ -380,7 +380,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj3 = z(i+3)
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER
+   !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER
    do j=1,$n
     S_inv(j,i  ) = S_inv(j,i  )*lambda - w(j)*zj
     S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1
@@ -435,15 +435,15 @@ BEGIN_TEMPLATE
 ! Version for mod(n,4) = 1
 subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m($n)  ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,$n)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   double precision               :: u($n), z($n), w($n), lambda, d_inv
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
   !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN
@@ -452,20 +452,20 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (LDS >= $n)
   integer                        :: i,j
   double precision :: zj, zj1, zj2, zj3
-  
+
   do i=1,$n
     u(i) = m(i) - S(i,l)
   enddo
 
   zj = 0.d0
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj)
+  !OMP$ SIMD REDUCTION(+:zj)
   do i=1,$n-1,4
     zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1)  &
-            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) 
+            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3)
   enddo
   zj = zj + S_inv($n,l)*u($n)
-  
+
   d_inv = 1.d0/d
   d = d+zj
   lambda = d*d_inv
@@ -473,7 +473,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
     d = 0.d0
     return
   endif
-  
+
   !DIR$ VECTOR ALIGNED
   do j=1,$n-1,4
    zj  = 0.d0
@@ -482,7 +482,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj3 = 0.d0
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
+   !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
    do i=1,$n-1
     zj  = zj  + S_inv(i,j  )*u(i)
     zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -498,19 +498,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj  = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER
+  !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER
   do i=1,$n-1
    zj = zj + S_inv(i,$n)*u(i)
   enddo
   z($n) = zj + S_inv($n,$n)*u($n)
 
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
+  !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
   do i=1,$n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
   enddo
-  
+
   do i=1,$n-1,4
    zj  = z(i  )
    zj1 = z(i+1)
@@ -518,14 +518,14 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj3 = z(i+3)
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER
+   !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3) NOVECREMAINDER
    do j=1,$n-1
     S_inv(j,i  ) = S_inv(j,i  )*lambda - w(j)*zj
     S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1
     S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2
     S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3
    enddo
-   S_inv($n,i  ) = S_inv($n,i  )*lambda - w($n)*zj 
+   S_inv($n,i  ) = S_inv($n,i  )*lambda - w($n)*zj
    S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*zj1
    S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*zj2
    S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*zj3
@@ -534,7 +534,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj = z($n)
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER
+  !OMP$ SIMD FIRSTPRIVATE(lambda,zj) NOVECREMAINDER
   do i=1,$n
    S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*zj
   enddo
@@ -588,15 +588,15 @@ BEGIN_TEMPLATE
 ! Version for mod(n,4) = 2
 subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m($n)  ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,$n)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   double precision               :: u($n), z($n), w($n), lambda, d_inv
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
   !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN
@@ -604,24 +604,24 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0)
   !DIR$ ASSUME (LDS >= $n)
   integer                        :: i,j
-  
+
   double precision :: zj, zj1, zj2, zj3
   !DIR$ NOPREFETCH
-  !DIR$ SIMD  NOVECREMAINDER
+  !OMP$ SIMD  NOVECREMAINDER
   do i=1,$n
     u(i) = m(i) - S(i,l)
   enddo
-  
+
   zj = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER
+  !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER
   do i=1,$n-2,4
     zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1)  &
-            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) 
+            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3)
   enddo
   i=$n-1
-  zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) 
+  zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1)
 
   d_inv = 1.d0/d
   d = d+zj
@@ -630,7 +630,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
     d = 0.d0
     return
   endif
-  
+
   !DIR$ VECTOR ALIGNED
   do j=1,$n-2,4
    zj  = 0.d0
@@ -638,7 +638,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj2 = 0.d0
    zj3 = 0.d0
    !DIR$ VECTOR ALIGNED
-   !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
+   !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3) NOVECREMAINDER
    do i=1,$n-2
     zj  = zj  + S_inv(i,j  )*u(i)
     zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -660,7 +660,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj1 = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER
+  !OMP$ SIMD REDUCTION(+:zj,zj1) NOVECREMAINDER
   do i=1,$n-2
    zj  = zj  + S_inv(i,j  )*u(i)
    zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -671,19 +671,19 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   z(j+1) = z(j+1) + S_inv($n,j+1)*u($n)
 
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
+  !OMP$ SIMD FIRSTPRIVATE(d_inv) NOVECREMAINDER
   do i=1,$n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
   enddo
-  
+
   do i=1,$n-2,4
    zj  = z(i)
    zj1 = z(i+1)
    zj2 = z(i+2)
    zj3 = z(i+3)
    !DIR$ VECTOR ALIGNED
-   !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)  NOVECREMAINDER
+   !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)  NOVECREMAINDER
    do j=1,$n-2
     S_inv(j,i  ) = S_inv(j,i  )*lambda -zj *w(j)
     S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j)
@@ -704,7 +704,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj = z(i)
   zj1= z(i+1)
   !DIR$ VECTOR ALIGNED
-  !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1)
+  !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1)
   do j=1,$n-2
    S_inv(j,i  ) = S_inv(j,i  )*lambda -zj*w(j)
    S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j)
@@ -761,15 +761,15 @@ BEGIN_TEMPLATE
 ! Version for mod(n,4) = 3
 subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS  ! Dimension of the vector
   real, intent(in)               :: m($n)  ! New vector
   integer, intent(in)            :: l      ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,$n)     ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix
   double precision,intent(inout) :: d            ! Det(S)
-  
+
   double precision               :: u($n), z($n), w($n), lambda, d_inv
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
   !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN
@@ -780,18 +780,18 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
 
   double precision :: zj, zj1, zj2, zj3
 
-  !DIR$ SIMD 
+  !OMP$ SIMD
   do i=1,$n
     u(i) = m(i) - S(i,l)
   enddo
-  
-  zj = 0.d0 
+
+  zj = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj) NOVECREMAINDER
+  !OMP$ SIMD REDUCTION(+:zj) NOVECREMAINDER
   do i=1,$n-3,4
     zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1)  &
-            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) 
+            + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3)
   enddo
   i=$n-2
   zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + S_inv(i+2,l)*u(i+2)
@@ -804,7 +804,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
     d = 0.d0
     return
   endif
-  
+
   !DIR$ VECTOR ALIGNED
   do j=1,$n-3,4
    zj  = 0.d0
@@ -812,7 +812,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj2 = 0.d0
    zj3 = 0.d0
    !DIR$ VECTOR ALIGNED
-   !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3)
+   !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3)
    do i=1,$n-3
     zj  = zj  + S_inv(i,j  )*u(i)
     zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -839,7 +839,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj2 = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zj,zj1,zj2)
+  !OMP$ SIMD REDUCTION(+:zj,zj1,zj2)
   do i=1,$n-3
    zj  = zj  + S_inv(i,j  )*u(i)
    zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -856,12 +856,12 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   z(j+2) = z(j+2) + S_inv($n,j+2)*u($n)
 
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv)
+  !OMP$ SIMD FIRSTPRIVATE(d_inv)
   do i=1,$n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
   enddo
-  
+
   do i=1,$n-3,4
    zj  = z(i)
    zj1 = z(i+1)
@@ -869,9 +869,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
    zj3 = z(i+3)
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)
+   !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)
    do j=1,$n-3
-    S_inv(j,i  ) = S_inv(j,i  )*lambda - w(j)*zj 
+    S_inv(j,i  ) = S_inv(j,i  )*lambda - w(j)*zj
     S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1
     S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2
     S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3
@@ -896,7 +896,7 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   zj2 = z(i+2)
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2)
+  !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2)
   do j=1,$n
    S_inv(j,i  ) = S_inv(j,i  )*lambda - w(j)*zj
    S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1
@@ -950,15 +950,15 @@ END_TEMPLATE
 
 subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
   implicit none
-  
+
   integer, intent(in)            :: n,LDS      ! Dimension of the vector
   real, intent(in)               :: m(LDS)     ! New vector
   integer, intent(in)            :: l          ! New position in S
-  
+
   real,intent(inout)             :: S(LDS,n)       ! Slater matrix
   double precision,intent(inout) :: S_inv(LDS,n)   ! Inverse Slater matrix
   double precision,intent(inout) :: d              ! Det(S)
-  
+
   double precision               :: lambda, d_inv
   double precision               :: u(3840), z(3840), w(3840)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u
@@ -976,19 +976,19 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
   do i=1,n
     u(i) = m(i) - S(i,l)
   enddo
-  
+
   zl = 0.d0
   !DIR$ VECTOR ALIGNED
   !DIR$ NOPREFETCH
-  !DIR$ SIMD REDUCTION(+:zl)
+  !OMP$ SIMD REDUCTION(+:zl)
   do i=1,n
    zl = zl + S_inv(i,l)*u(i)
   enddo
 
-  d_inv = 1.d0/d 
+  d_inv = 1.d0/d
   d = d+zl
   lambda = d*d_inv
-  
+
   if ( dabs(lambda) < 1.d-3 ) then
     d = 0.d0
     return
@@ -1006,7 +1006,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
    zj3 = 0.d0
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD REDUCTION(+:zj,zj1,zj2,zj3)
+   !OMP$ SIMD REDUCTION(+:zj,zj1,zj2,zj3)
    do i=1,n
     zj  = zj  + S_inv(i,j  )*u(i)
     zj1 = zj1 + S_inv(i,j+1)*u(i)
@@ -1023,7 +1023,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
    zj = 0.d0
    !DIR$ VECTOR ALIGNED
    !DIR$ NOPREFETCH
-   !DIR$ SIMD REDUCTION(+:zj)
+   !OMP$ SIMD REDUCTION(+:zj)
    do i=1,n
     zj = zj + S_inv(i,j)*u(i)
    enddo
@@ -1031,14 +1031,14 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
   enddo
 
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv)
+  !OMP$ SIMD FIRSTPRIVATE(d_inv)
   do i=1,n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
   enddo
-  
+
   !DIR$ NOPREFETCH
-  !DIR$ SIMD FIRSTPRIVATE(d_inv)
+  !OMP$ SIMD FIRSTPRIVATE(d_inv)
   do i=1,n
     w(i) = S_inv(i,l)*d_inv
     S(i,l) = m(i)
@@ -1051,7 +1051,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
     zj3 = z(i+3)
     !DIR$ VECTOR ALIGNED
     !DIR$ NOPREFETCH
-    !DIR$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)
+    !OMP$ SIMD FIRSTPRIVATE(lambda,zj,zj1,zj2,zj3)
     do j=1,n
      S_inv(j,i  ) = S_inv(j,i  )*lambda -zj *w(j)
      S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j)
@@ -1064,7 +1064,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
     zj = z(i)
     !DIR$ VECTOR ALIGNED
     !DIR$ NOPREFETCH
-    !DIR$ SIMD FIRSTPRIVATE(lambda,zj)
+    !OMP$ SIMD FIRSTPRIVATE(lambda,zj)
     do j=1,n
       S_inv(j,i) = S_inv(j,i)*lambda -zj*w(j)
     enddo
@@ -1146,9 +1146,9 @@ END_PROVIDER
  BEGIN_PROVIDER [ double precision, det_alpha_value_curr ]
 &BEGIN_PROVIDER [ real, slater_matrix_alpha, (elec_alpha_num_8,elec_alpha_num) ]
 &BEGIN_PROVIDER [ double precision, slater_matrix_alpha_inv_det, (elec_alpha_num_8,elec_alpha_num) ]
-  
+
   implicit none
-  
+
   BEGIN_DOC
   ! det_alpha_value_curr : Value of the current alpha determinant
   !
@@ -1159,7 +1159,7 @@ END_PROVIDER
   !
   ! slater_matrix_alpha_inv_det: Inverse of the alpha Slater matrix x determinant
   END_DOC
-  
+
   double precision               :: ddet
   integer                        :: i,j,k,imo,l
   integer                        :: to_do(elec_alpha_num), n_to_do_old, n_to_do
@@ -1172,7 +1172,7 @@ END_PROVIDER
     !DIR$ VECTOR ALIGNED
     slater_matrix_alpha_inv_det = 0.d0
   endif
-  
+
   PROVIDE mo_value
   if (det_i /= det_alpha_order(1) ) then
 !  if (det_i == -1 ) then
@@ -1187,7 +1187,7 @@ END_PROVIDER
     enddo
 
     ddet = 0.d0
-    
+
     if (n_to_do < shiftl(elec_alpha_num,1)) then
 
       do while ( n_to_do > 0 )
@@ -1218,13 +1218,13 @@ END_PROVIDER
       enddo
 
     endif
-    
+
   else
 
     ddet = 0.d0
 
   endif
-    
+
   ! Avoid NaN
   if (ddet /= 0.d0) then
     continue
@@ -1246,11 +1246,11 @@ END_PROVIDER
       enddo
     enddo
     call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet)
-    
+
   endif
   ASSERT (ddet /= 0.d0)
-  
-  det_alpha_value_curr = ddet 
+
+  det_alpha_value_curr = ddet
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
@@ -1266,24 +1266,24 @@ END_PROVIDER
   !
   !  slater_matrix_beta_inv_det : Inverse of the beta Slater matrix x determinant
   END_DOC
-  
+
   double precision               :: ddet
   integer                        :: i,j,k,imo,l
   integer                        :: to_do(elec_alpha_num-mo_closed_num), n_to_do_old, n_to_do
-  
+
   integer, save                  :: ifirst
   if (elec_beta_num == 0) then
     det_beta_value_curr = 0.d0
     return
   endif
-  
+
   if (ifirst == 0) then
     ifirst = 1
     slater_matrix_beta = 0.
     slater_matrix_beta_inv_det = 0.d0
   endif
   PROVIDE mo_value
-  
+
   if (det_j /= det_beta_order(1)) then
 !  if (det_j == -1) then
 
@@ -1295,7 +1295,7 @@ END_PROVIDER
           to_do(n_to_do) = k
       endif
     enddo
-    
+
     ddet = 0.d0
     if (n_to_do < shiftl(elec_beta_num,1)) then
 
@@ -1327,13 +1327,13 @@ END_PROVIDER
       enddo
 
     endif
-    
+
   else
-    
+
     ddet = 0.d0
 
   endif
-  
+
   ! Avoid NaN
   if (ddet /= 0.d0) then
     continue
@@ -1357,9 +1357,9 @@ END_PROVIDER
     call invert(slater_matrix_beta_inv_det,elec_beta_num_8,elec_beta_num,ddet)
   endif
   ASSERT (ddet /= 0.d0)
-  
-  det_beta_value_curr = ddet 
-  
+
+  det_beta_value_curr = ddet
+
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer, det_alpha_num_pseudo ]
@@ -1381,15 +1381,15 @@ END_PROVIDER
  BEGIN_PROVIDER [ double precision , det_alpha_value,  (det_alpha_num_8) ]
 &BEGIN_PROVIDER [ double precision , det_alpha_grad_lapl, (4,elec_alpha_num,det_alpha_num) ]
 &BEGIN_PROVIDER [ double precision , det_alpha_pseudo, (elec_alpha_num_8,det_alpha_num_pseudo) ]
-  
+
   implicit none
-  
+
   BEGIN_DOC
   ! Values of the alpha determinants
   ! Gradients of the alpha determinants
   ! Laplacians of the alpha determinants
   END_DOC
-  
+
   integer                        :: j,i,k
   integer, save                  :: ifirst = 0
   if (ifirst == 0) then
@@ -1398,73 +1398,73 @@ END_PROVIDER
     det_alpha_grad_lapl = 0.d0
     det_alpha_pseudo = 0.d0
   endif
-  
-  
+
+
   do j=1,det_alpha_num
-    
+
     det_i = det_alpha_order(j)
     if (j > 1) then
       TOUCH det_i
     endif
-    
+
     det_alpha_value(det_i) = det_alpha_value_curr
     det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:)
     if (do_pseudo) then
-      det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) 
+      det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:)
     endif
-   
+
   enddo
-  
+
   det_i = det_alpha_order(1)
   SOFT_TOUCH det_i
-  
+
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, det_beta_value,  (det_beta_num_8) ]
 &BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_alpha_num+1:elec_num,det_beta_num) ]
 &BEGIN_PROVIDER [ double precision, det_beta_pseudo, (elec_alpha_num+1:elec_num,det_beta_num_pseudo) ]
-  
-  
+
+
   implicit none
-  
+
   BEGIN_DOC
   ! Values of the beta determinants
   ! Gradients of the beta determinants
   ! Laplacians of the beta determinants
   END_DOC
-  
+
   integer                        :: j,i,k
   integer, save                  :: ifirst = 0
   if (elec_beta_num == 0) then
     det_beta_value = 1.d0
     return
   endif
-    
+
   if (ifirst == 0) then
     ifirst = 1
     det_beta_value  = 0.d0
     det_beta_grad_lapl = 0.d0
     det_beta_pseudo = 0.d0
   endif
-  
+
   do j=1,det_beta_num
-    
+
     det_j = det_beta_order(j)
     if (j > 1) then
       TOUCH det_j
     endif
-    
+
     det_beta_value(det_j)  = det_beta_value_curr
     det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:)
     if (do_pseudo) then
-      det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) 
+      det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:)
     endif
 
   enddo
-  
+
   det_j = det_beta_order(1)
   SOFT_TOUCH det_j
-  
+
 END_PROVIDER
 
 
@@ -1472,34 +1472,34 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, psidet_inv ]
 &BEGIN_PROVIDER [ double precision, psidet_grad_lapl, (4,elec_num_8) ]
 &BEGIN_PROVIDER [ double precision, pseudo_non_local, (elec_num) ]
-  
+
   implicit none
   BEGIN_DOC
   ! Value of the determinantal part of the wave function
-  
+
   ! Gradient of the determinantal part of the wave function
-  
+
   ! Laplacian of determinantal part of the wave function
 
   ! Non-local component of the pseudopotentials
 
   ! Regularized 1/psi = 1/(psi + 1/psi)
   END_DOC
-  
+
   integer, save                  :: ifirst=0
   if (ifirst == 0) then
     ifirst = 1
     psidet_grad_lapl = 0.d0
   endif
-  
+
   double precision :: CDb(det_alpha_num_8)
   double precision :: CDb_i
   double precision :: DaC(det_beta_num_8)
   !DIR$ ATTRIBUTES ALIGN : 32 :: DaC,CDb
 
   ! C x D_beta
-  ! D_alpha^t x C 
-  ! D_alpha^t x (C x D_beta) 
+  ! D_alpha^t x C
+  ! D_alpha^t x (C x D_beta)
 
   integer :: i,j,k, l
   integer :: i1,i2,i3,i4,det_num4
@@ -1576,7 +1576,7 @@ END_PROVIDER
     psidet_value = psidet_value + det_beta_value(j) * DaC(j)
   enddo
 
-  
+
   if (psidet_value == 0.d0) then
     call abrt(irp_here,'Determinantal component of the wave function is zero.')
   endif
@@ -1595,7 +1595,7 @@ END_PROVIDER
         size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2),         &
         DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1)
   endif
-  
+
   if (do_pseudo) then
     call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv,          &
         det_alpha_pseudo, size(det_alpha_pseudo,1),                  &
@@ -1669,7 +1669,7 @@ BEGIN_PROVIDER  [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
   BEGIN_DOC
   ! Gradient of the current alpha determinant
   END_DOC
-  
+
   integer                        :: i, j, k
   !DIR$ VECTOR ALIGNED
   do i=1,elec_alpha_num
@@ -1678,7 +1678,7 @@ BEGIN_PROVIDER  [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
     det_alpha_grad_lapl_curr(3,i) = 0.d0
     det_alpha_grad_lapl_curr(4,i) = 0.d0
   enddo
-  
+
   integer :: imo, imo2
 
 ! -------
@@ -1688,7 +1688,7 @@ BEGIN_PROVIDER  [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
 !      imo  = mo_list_alpha_curr(j)
 !      do i=1,elec_alpha_num
 !        do k=1,4
-!          det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j) 
+!          det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j)
 !        enddo
 !      enddo
 !    enddo
@@ -1738,13 +1738,13 @@ BEGIN_PROVIDER  [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
     do i=1,elec_alpha_num
         !DIR$ VECTOR ALIGNED
       do k=1,4
-        det_alpha_grad_lapl_curr(k,i  ) = det_alpha_grad_lapl_curr(k,i  ) + mo_grad_lapl_alpha(k,i  ,imo)*slater_matrix_alpha_inv_det(i  ,j) 
+        det_alpha_grad_lapl_curr(k,i  ) = det_alpha_grad_lapl_curr(k,i  ) + mo_grad_lapl_alpha(k,i  ,imo)*slater_matrix_alpha_inv_det(i  ,j)
       enddo
     enddo
 
   endif
-    
-  
+
+
 END_PROVIDER
 
 
@@ -1753,7 +1753,7 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
   BEGIN_DOC
   ! Gradient and Laplacian of the current beta determinant
   END_DOC
-  
+
   integer                        :: i, j, k, l
 
   !DIR$ VECTOR ALIGNED
@@ -1763,7 +1763,7 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
     det_beta_grad_lapl_curr(3,i) = 0.d0
     det_beta_grad_lapl_curr(4,i) = 0.d0
   enddo
-  
+
   integer                          :: imo, imo2
 
 ! -------
@@ -1779,7 +1779,7 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
 !    enddo
 !  enddo
 !
-! -------- 
+! --------
 
   if (iand(elec_beta_num,1) == 0) then
 
@@ -1793,10 +1793,10 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
         do k=1,4
           det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +&
               mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j  ) + &
-              mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) 
+              mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1)
           det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +&
               mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j  ) + &
-              mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) 
+              mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1)
         enddo
       enddo
     enddo
@@ -1813,10 +1813,10 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
         do k=1,4
           det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +&
               mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j  ) + &
-              mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) 
+              mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1)
           det_beta_grad_lapl_curr(k,i+1) = det_beta_grad_lapl_curr(k,i+1) +&
               mo_grad_lapl_beta(k,i+1,imo )*slater_matrix_beta_inv_det(l+1,j  ) + &
-              mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1) 
+              mo_grad_lapl_beta(k,i+1,imo2)*slater_matrix_beta_inv_det(l+1,j+1)
         enddo
       enddo
       i = elec_num
@@ -1825,7 +1825,7 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
       do k=1,4
         det_beta_grad_lapl_curr(k,i) = det_beta_grad_lapl_curr(k,i) +&
             mo_grad_lapl_beta(k,i,imo )*slater_matrix_beta_inv_det(l,j  ) + &
-            mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1) 
+            mo_grad_lapl_beta(k,i,imo2)*slater_matrix_beta_inv_det(l,j+1)
       enddo
     enddo
 
@@ -1840,7 +1840,7 @@ BEGIN_PROVIDER  [ double precision, det_beta_grad_lapl_curr, (4,elec_alpha_num+1
       enddo
     enddo
 
-  endif      
+  endif
 
 END_PROVIDER
 
diff --git a/src/mo.irp.f b/src/mo.irp.f
index 38228b2..124890d 100644
--- a/src/mo.irp.f
+++ b/src/mo.irp.f
@@ -5,12 +5,12 @@
 ! Number of Molecular orbitals
   END_DOC
   integer, external              :: mod_align
-  
+
   mo_num = maxval(present_mos)
   call iinfo(irp_here,'mo_num',mo_num)
-  
+
   mo_num_8 = mod_align(mo_num)
-  
+
 END_PROVIDER
 
 
@@ -22,7 +22,7 @@ BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ]
   integer                        :: i, j
   real,allocatable               :: buffer(:,:)
   allocate (buffer(ao_num,mo_tot_num))
-  
+
   buffer = 0.
   call get_mo_basis_mo_coef(buffer)
   do i=1,mo_tot_num
@@ -35,7 +35,7 @@ BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ]
     call set_order(mo_coef_input(1,i),ao_nucl_sort_idx,ao_num)
   enddo
   deallocate(buffer)
-  
+
 END_PROVIDER
 
 
@@ -56,7 +56,7 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ]
 ! Molecular orbital coefficients
   END_DOC
   integer                        :: i, j
-  
+
   do j=1,mo_num
     do i=1,ao_num_8
       mo_coef(i,j) = mo_coef_input(i,j)
@@ -68,10 +68,10 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ]
       mo_coef(i,j) = 0.
     enddo
   enddo
-  
+
   ! Input MOs are not needed any more
   FREE mo_coef_input
-  
+
   real                           :: f
   f = 1./mo_scale
   do j=1,mo_num
@@ -80,7 +80,7 @@ BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ]
       mo_coef(i,j) *= f
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -90,7 +90,7 @@ BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num_8,ao_num_8) ]
 ! Transpose of the Molecular orbital coefficients
   END_DOC
   call transpose(mo_coef,ao_num_8,mo_coef_transp,mo_num_8,ao_num_8,mo_num_8)
-  
+
 END_PROVIDER
 
 
@@ -102,7 +102,7 @@ END_PROVIDER
 ! orbital coefficients
   END_DOC
   integer                        :: i, j
-  
+
   integer                        :: idx
   mo_coef_transp_sparsity = 0.
   do j=1,ao_num
@@ -117,7 +117,7 @@ END_PROVIDER
     mo_coef_transp_sparsity += float(idx)
   enddo
   mo_coef_transp_sparsity *= 1./(mo_num*ao_num)
-  
+
 END_PROVIDER
 
 
@@ -133,7 +133,7 @@ BEGIN_PROVIDER [ real, mo_coef_transp_present, (num_present_mos_8,ao_num_8) ]
       mo_coef_transp_present(j,i) = mo_coef_transp(present_mos(j),i)
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -143,24 +143,24 @@ END_PROVIDER
 &BEGIN_PROVIDER [ real, mo_grad_transp_z, (mo_num_8,elec_num) ]
 &BEGIN_PROVIDER [ real, mo_lapl_transp, (mo_num_8,elec_num) ]
   implicit none
-  
+
   BEGIN_DOC
 ! Values, gradients, laplacians of the molecular orbitals
 !
 ! Arrays are padded for efficiency
   END_DOC
-  
+
   integer                        :: i, j, k, l, m
-  
+
   PROVIDE primitives_reduced
-  
+
   if (do_nucl_fitcusp) then
     PROVIDE nucl_fitcusp_param
     PROVIDE nucl_elec_dist_vec
     PROVIDE nucl_elec_dist_inv
   endif
-  
-  
+
+
   do i=1,elec_num
     if (i>1) then
       ao_elec = i
@@ -180,7 +180,7 @@ END_PROVIDER
           mo_grad_transp_z(1,i),                                     &
           mo_lapl_transp(1,i),                                       &
           ao_num)
-      
+
     else
       call sparse_full_mv(mo_coef_transp_present,num_present_mos_8,  &
           ao_value_block(1),ao_num_8,                                &
@@ -195,7 +195,7 @@ END_PROVIDER
           mo_grad_transp_z(1,i),                                     &
           mo_lapl_transp(1,i),                                       &
           ao_num)
-      
+
       do j=num_present_mos,1,-1
         mo_value_transp (present_mos(j),i) = mo_value_transp (j,i)
         mo_grad_transp_x(present_mos(j),i) = mo_grad_transp_x(j,i)
@@ -207,10 +207,10 @@ END_PROVIDER
         endif
       enddo
     endif
-    
+
     if (do_nucl_fitcusp) then
       real                           :: r, r2, r_inv, d, expzr, Z, Z2, a, b, c, phi, rx, ry, rz
-      
+
       do k=1,nucl_num
         r = nucl_elec_dist(k,i)
         if (r > nucl_fitcusp_radius(k)) then
@@ -235,14 +235,14 @@ END_PROVIDER
         enddo
         exit
       enddo ! k
-      
+
     endif
-    
+
   enddo  ! i
   ao_elec = 1
   SOFT_TOUCH ao_elec
-  
-  
+
+
   if (do_prepare) then
     real                           :: lambda, t
     ! Scale off-diagonal elements
@@ -268,7 +268,7 @@ END_PROVIDER
       enddo
     enddo
   endif
-  
+
   do i=1,mo_num
     do j=1,elec_num
         mo_value_transp(i,j)  *= mo_cusp_rescale(i)
@@ -285,10 +285,10 @@ BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ]
   BEGIN_DOC
 ! Values of the molecular orbitals
   END_DOC
-  
+
   integer                        :: i,j
   integer, save                  :: ifirst = 0
-  
+
   if (ifirst == 0) then
     ifirst = 1
     PROVIDE primitives_reduced
@@ -296,7 +296,7 @@ BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ]
     mo_value = 0.
   endif
   call transpose(mo_value_transp(1,1),mo_num_8,mo_value,elec_num_8,mo_num,elec_num)
-  
+
 END_PROVIDER
 
 
@@ -304,14 +304,14 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, mo_grad_y, (elec_num_8,mo_num) ]
 &BEGIN_PROVIDER [ double precision, mo_grad_z, (elec_num_8,mo_num) ]
   implicit none
-  
+
   BEGIN_DOC
 ! Gradients of the molecular orbitals
   END_DOC
-  
+
   integer                        :: i,j
   integer, save                  :: ifirst = 0
-  
+
   if (ifirst == 0) then
     !DIR$ VECTOR ALIGNED
     mo_grad_x = 0.d0
@@ -327,7 +327,7 @@ END_PROVIDER
   call transpose_to_dp(mo_grad_transp_z(1,1),mo_num_8,mo_grad_z(1,1),elec_num_8,mo_num,elec_num)
   call transpose_to_dp(mo_grad_transp_x(1,1),mo_num_8,mo_grad_x(1,1),elec_num_8,mo_num,elec_num)
 
-  
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ]
@@ -335,7 +335,7 @@ BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ]
   BEGIN_DOC
 ! Laplacians of the molecular orbitals
   END_DOC
-  
+
   integer                        :: i,j
   integer, save                  :: ifirst = 0
   if (ifirst == 0) then
@@ -345,7 +345,7 @@ BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ]
     mo_lapl = 0.d0
   endif
   call transpose_to_dp(mo_lapl_transp(1,1),mo_num_8,mo_lapl,elec_num_8,mo_num,elec_num)
-  
+
 END_PROVIDER
 
 
@@ -358,21 +358,21 @@ END_PROVIDER
  integer :: i,j
  do j=1,mo_num
     do i=1,elec_alpha_num
-      mo_grad_lapl_alpha(1,i,j) = mo_grad_transp_x(j,i) 
-      mo_grad_lapl_alpha(2,i,j) = mo_grad_transp_y(j,i) 
-      mo_grad_lapl_alpha(3,i,j) = mo_grad_transp_z(j,i) 
-      mo_grad_lapl_alpha(4,i,j) = mo_lapl_transp  (j,i) 
+      mo_grad_lapl_alpha(1,i,j) = mo_grad_transp_x(j,i)
+      mo_grad_lapl_alpha(2,i,j) = mo_grad_transp_y(j,i)
+      mo_grad_lapl_alpha(3,i,j) = mo_grad_transp_z(j,i)
+      mo_grad_lapl_alpha(4,i,j) = mo_lapl_transp  (j,i)
     enddo
   enddo
   do j=1,mo_num
     do i=elec_alpha_num+1,elec_num
-      mo_grad_lapl_beta(1,i,j) = mo_grad_transp_x(j,i) 
-      mo_grad_lapl_beta(2,i,j) = mo_grad_transp_y(j,i) 
-      mo_grad_lapl_beta(3,i,j) = mo_grad_transp_z(j,i) 
-      mo_grad_lapl_beta(4,i,j) = mo_lapl_transp  (j,i) 
+      mo_grad_lapl_beta(1,i,j) = mo_grad_transp_x(j,i)
+      mo_grad_lapl_beta(2,i,j) = mo_grad_transp_y(j,i)
+      mo_grad_lapl_beta(3,i,j) = mo_grad_transp_z(j,i)
+      mo_grad_lapl_beta(4,i,j) = mo_lapl_transp  (j,i)
   enddo
  enddo
- 
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_grad_lapl_transp, (4,mo_num,elec_num) ]
@@ -389,7 +389,7 @@ BEGIN_PROVIDER [ double precision, mo_grad_lapl_transp, (4,mo_num,elec_num) ]
      mo_grad_lapl_transp(4,j,i) = mo_lapl_transp  (j,i)
    enddo
  enddo
- 
+
 END_PROVIDER
 
 
@@ -404,18 +404,18 @@ END_PROVIDER
 
 
 BEGIN_PROVIDER [ integer, mo_tot_num ]
-  
+
   BEGIN_DOC
 ! Total number of MOs in the EZFIO file
   END_DOC
-  
+
   mo_tot_num = -1
   call get_mo_basis_mo_num(mo_tot_num)
   if (mo_tot_num <= 0) then
     call abrt(irp_here,'Total number of MOs can''t be <0')
   endif
   call iinfo(irp_here,'mo_tot_num',mo_tot_num)
-  
+
 END_PROVIDER
 
 
@@ -432,16 +432,16 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
   END_DOC
   integer                        :: i, j, k, l
   real                           :: ao_value_at_nucl_no_S(ao_num)
-  
+
   PROVIDE mo_fitcusp_normalization_before
   do k=1,nucl_num
     point(1) = nucl_coord(k,1)
     point(2) = nucl_coord(k,2)
     point(3) = nucl_coord(k,3)
     TOUCH point
-    
+
     PROVIDE ao_value_p
-    
+
     do i=1,ao_num
       if (ao_nucl(i) /= k) then
         ao_value_at_nucl_no_S(i) = ao_value_p(i)
@@ -449,7 +449,7 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
         ao_value_at_nucl_no_S(i) = 0.
       endif
     enddo
-    
+
     integer :: jj
     do jj=1,num_present_mos
       j = present_mos(jj)
@@ -459,11 +459,11 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
         mo_value_at_nucl(j,k) = mo_value_at_nucl(j,k) + mo_coef(i,j)*ao_value_at_nucl_no_S(i)
       enddo
     enddo
-    
+
   enddo
   FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
   SOFT_TOUCH point
-  
+
 END_PROVIDER
 
 
@@ -474,15 +474,15 @@ END_PROVIDER
   BEGIN_DOC
 ! Values of the atomic orbitals without S components on atoms
   END_DOC
-  
+
   integer                        :: i, j, k
-  
+
   do k=1,nucl_num
     point(1) = nucl_coord(k,1)
     point(2) = nucl_coord(k,2)
     point(3) = nucl_coord(k,3)+ nucl_fitcusp_radius(k)
     TOUCH point
-    
+
     do j=1,ao_num
       ao_value_at_fitcusp_radius(j,k) = ao_value_p(j)
       ao_grad_at_fitcusp_radius(j,k) = ao_grad_p(j,3)
@@ -496,7 +496,7 @@ END_PROVIDER
   enddo
   FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
   SOFT_TOUCH point
-  
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_before, (mo_tot_num) ]
@@ -560,7 +560,7 @@ BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num)
           t = 0.d0
           do j=1,ao_num
             if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then
-              t = t + mo_coef(j,i) * ao_value_p(j) 
+              t = t + mo_coef(j,i) * ao_value_p(j)
             endif
           enddo
           t = t + nucl_fitcusp_param(1,i,k) +                        &
@@ -601,7 +601,7 @@ END_PROVIDER
 ! Values of the molecular orbitals without S components on atoms
   END_DOC
   integer                        :: i, j, k, l
-  
+
   do k=1,nucl_num
     do j=1,mo_num
       mo_value_at_fitcusp_radius(j,k) = 0.d0
@@ -615,7 +615,7 @@ END_PROVIDER
       enddo
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -626,7 +626,7 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ]
   END_DOC
   integer                        :: i,k, niter
   character*(80)                 :: message
-  
+
   nucl_fitcusp_param = 0.d0
   do k=1,nucl_num
     double precision               :: r, Z
@@ -642,14 +642,14 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ]
       grad_phi = mo_grad_at_fitcusp_radius(i,k)
       phi = mo_value_at_fitcusp_radius(i,k)
       eta = mo_value_at_nucl(i,k)
-      
+
       nucl_fitcusp_param(1,i,k) = -(R*(2.d0*eta*Z-6.d0*grad_phi)+lap_phi*R*R+6.d0*phi)/(2.d0*R*Z-6.d0)
       nucl_fitcusp_param(2,i,k) = (lap_phi*R*R*Z-6.d0*grad_phi*R*Z+6.d0*phi*Z+6.d0*eta*Z)/(2.d0*R*Z-6.d0)
       nucl_fitcusp_param(3,i,k) = -(R*(-5.d0*grad_phi*Z-1.5d0*lap_phi)+lap_phi*R*R*Z+3.d0*phi*Z+&
           3.d0*eta*Z+6.d0*grad_phi)/(R*R*Z-3.d0*R)
       nucl_fitcusp_param(4,i,k) = (R*(-2.d0*grad_phi*Z-lap_phi)+0.5d0*lap_phi*R*R*Z+phi*Z+&
           eta*Z+3.d0*grad_phi)/(R*R*R*Z-3.d0*R*R)
-      
+
     enddo
   enddo
 END_PROVIDER
@@ -689,7 +689,7 @@ subroutine sparse_full_mv(A,LDA,                                     &
   !DIR$ ASSUME_ALIGNED C3 : $IRP_ALIGN
   !DIR$ ASSUME_ALIGNED C4 : $IRP_ALIGN
   !DIR$ ASSUME_ALIGNED C5 : $IRP_ALIGN
-  
+
   integer                        :: kao, kmax, kmax2, kmax3
   integer                        :: i,j,k
   integer                        :: k_vec(8)
@@ -698,9 +698,9 @@ subroutine sparse_full_mv(A,LDA,                                     &
   real                           :: d21, d22, d23, d24, d25
   real                           :: d31, d32, d33, d34, d35
   real                           :: d41, d42, d43, d44, d45
-  
+
   ! LDC and LDA have to be factors of simd_sp
-  
+
 !  IRP_IF NO_PREFETCH
 !  IRP_ELSE
 !    call MM_PREFETCH (A(1,indices(1)),3)
@@ -709,7 +709,7 @@ subroutine sparse_full_mv(A,LDA,                                     &
 !    call MM_PREFETCH (A(1,indices(4)),3)
 !  IRP_ENDIF
 
-  !DIR$ SIMD
+  !OMP$ SIMD
   do j=1,LDC
     C1(j) = 0.
     C2(j) = 0.
@@ -727,35 +727,35 @@ subroutine sparse_full_mv(A,LDA,                                     &
     k_vec(2) = indices(kao+1)
     k_vec(3) = indices(kao+2)
     k_vec(4) = indices(kao+3)
-    
+
     d11 = B1(kao  )
     d21 = B1(kao+1)
     d31 = B1(kao+2)
     d41 = B1(kao+3)
-    
+
     d12 = B2(kao  )
     d22 = B2(kao+1)
     d32 = B2(kao+2)
     d42 = B2(kao+3)
-    
+
     d13 = B3(kao  )
     d23 = B3(kao+1)
     d33 = B3(kao+2)
     d43 = B3(kao+3)
-    
+
     d14 = B4(kao  )
     d24 = B4(kao+1)
     d34 = B4(kao+2)
     d44 = B4(kao+3)
-    
+
     d15 = B5(kao  )
     d25 = B5(kao+1)
     d35 = B5(kao+2)
     d45 = B5(kao+3)
-    
+
     do k=0,LDA-1,$IRP_ALIGN/4
       !DIR$ VECTOR ALIGNED
-      !DIR$ SIMD FIRSTPRIVATE(d11,d21,d31,d41)
+      !OMP$ SIMD FIRSTPRIVATE(d11,d21,d31,d41)
       do j=1,$IRP_ALIGN/4
 !        IRP_IF NO_PREFETCH
 !        IRP_ELSE
@@ -767,18 +767,18 @@ subroutine sparse_full_mv(A,LDA,                                     &
         C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 + A(j+k,k_vec(2))*d21&
             +   A(j+k,k_vec(3))*d31 + A(j+k,k_vec(4))*d41
       enddo
-    
+
       !DIR$ VECTOR ALIGNED
-      !DIR$ SIMD FIRSTPRIVATE(d12,d22,d32,d42,d13,d23,d33,d43)
+      !OMP$ SIMD FIRSTPRIVATE(d12,d22,d32,d42,d13,d23,d33,d43)
       do j=1,$IRP_ALIGN/4
         C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12 + A(j+k,k_vec(2))*d22&
             +   A(j+k,k_vec(3))*d32 + A(j+k,k_vec(4))*d42
         C3(j+k) = C3(j+k) + A(j+k,k_vec(1))*d13 + A(j+k,k_vec(2))*d23&
             +   A(j+k,k_vec(3))*d33 + A(j+k,k_vec(4))*d43
       enddo
-    
+
       !DIR$ VECTOR ALIGNED
-      !DIR$ SIMD FIRSTPRIVATE(d14,d24,d34,d44,d15,d25,d35,d45)
+      !OMP$ SIMD FIRSTPRIVATE(d14,d24,d34,d44,d15,d25,d35,d45)
       do j=1,$IRP_ALIGN/4
         C4(j+k) = C4(j+k) + A(j+k,k_vec(1))*d14 + A(j+k,k_vec(2))*d24&
             +   A(j+k,k_vec(3))*d34 + A(j+k,k_vec(4))*d44
@@ -786,9 +786,9 @@ subroutine sparse_full_mv(A,LDA,                                     &
             +   A(j+k,k_vec(3))*d35 + A(j+k,k_vec(4))*d45
       enddo
     enddo
-    
+
   enddo
-  
+
   do kao = kmax2+1, kmax3
     k_vec(1) = indices(kao)
     d11 = B1(kao)
@@ -799,7 +799,7 @@ subroutine sparse_full_mv(A,LDA,                                     &
     !DIR$ VECTOR ALIGNED
     do k=0,LDA-1,$IRP_ALIGN/4
       !DIR$ VECTOR ALIGNED
-      !DIR$ SIMD FIRSTPRIVATE(d11,d12,d13,d14,d15)
+      !OMP$ SIMD FIRSTPRIVATE(d11,d12,d13,d14,d15)
       do j=1,$IRP_ALIGN/4
         C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11
         C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12
@@ -809,7 +809,7 @@ subroutine sparse_full_mv(A,LDA,                                     &
       enddo
     enddo
   enddo
-  
+
 end