From 0cc458a23f07bdd94ac9faefa23dc0ac2cd3a98f Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Tue, 19 Jun 2018 11:14:44 +0200
Subject: [PATCH] Occupation numbers in EZFIO

---
 lib64                                 |  1 -
 ocaml/qptypes_generator.ml            |  4 +-
 src/Determinants/density_matrix.irp.f |  3 +-
 src/MO_Basis/utils.irp.f              | 59 +++++++++++++++++++++++++++
 4 files changed, 63 insertions(+), 4 deletions(-)
 delete mode 120000 lib64

diff --git a/lib64 b/lib64
deleted file mode 120000
index 7951405f..00000000
--- a/lib64
+++ /dev/null
@@ -1 +0,0 @@
-lib
\ No newline at end of file
diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml
index d76120b8..3d5950a5 100644
--- a/ocaml/qptypes_generator.ml
+++ b/ocaml/qptypes_generator.ml
@@ -84,8 +84,8 @@ let input_data = "
 * MO_coef : float
 
 * MO_occ : float
-  if (x < 0.) || (x > 2.) then
-    raise (Invalid_argument (Printf.sprintf \"MO_occ : (0. <= x <= 2.) : x=%f\" x)); 
+  if x < 0. then 0.  else 
+  if x > 2. then 2.  else
 
 * AO_coef : float
 
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index 0bd0a187..ea32ece4 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -347,7 +347,8 @@ subroutine set_natural_mos
  double precision, allocatable :: tmp(:,:)
 
  label = "Natural"
- call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label)
+ call mo_as_svd_vectors_of_mo_matrix_eig(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,mo_occ,label)
+ soft_touch mo_occ
 
 end
 subroutine save_natural_mos
diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f
index b8eed315..3882aae4 100644
--- a/src/MO_Basis/utils.irp.f
+++ b/src/MO_Basis/utils.irp.f
@@ -157,6 +157,65 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   mo_label = label
 end
 
+subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
+  implicit none
+  integer,intent(in)             :: lda,m,n
+  character*(64), intent(in)     :: label
+  double precision, intent(in)   :: matrix(lda,n)
+  double precision, intent(out)  :: eig(m)
+  
+  integer :: i,j
+  double precision :: accu
+  double precision, allocatable  :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:)
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
+  
+  call write_time(6)
+  if (m /= mo_tot_num) then
+    print *, irp_here, ': Error : m/= mo_tot_num'
+    stop 1
+  endif
+
+  allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num,m),D(m),Vt(lda,n))
+
+  do j=1,n
+    do i=1,m
+      A(i,j) = matrix(i,j)
+    enddo
+  enddo
+  mo_coef_new = mo_coef
+  
+  call svd(A,lda,U,lda,D,Vt,lda,m,n)
+
+  write (6,'(A)') 'MOs are now **'//trim(label)//'**'
+  write (6,'(A)')  ''
+  write (6,'(A)') 'Eigenvalues'
+  write (6,'(A)')  '-----------'
+  write (6,'(A)') ''
+  write (6,'(A)')  '======== ================ ================'
+  write (6,'(A)')  '   MO       Eigenvalue       Cumulative   '
+  write (6,'(A)')  '======== ================ ================'
+
+  accu = 0.d0
+  do i=1,m
+    accu = accu + D(i)
+    write (6,'(I8,1X,F16.10,1X,F16.10)')  i,D(i), accu
+  enddo
+  write (6,'(A)')  '======== ================ ================'
+  write (6,'(A)')  ''
+  
+  call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
+
+  do i=1,m
+    eig(i) = D(i)
+  enddo
+
+  deallocate(A,mo_coef_new,U,Vt,D)
+  call write_time(6)
+  
+  mo_label = label
+
+end
+
 subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
   implicit none
   integer,intent(in)             :: n,m