Merge lpqsv26:~/quantum_package

lib64

@@ -1 +0,0 @@
-lib

ocaml/qptypes_generator.ml

@@ -84,8 +84,8 @@ let input_data = "
 * MO_coef : float
 
 * MO_occ : float
-  if (x < 0.) || (x > 2.) then
+  if x < 0. then 0.  else 
-    raise (Invalid_argument (Printf.sprintf \"MO_occ : (0. <= x <= 2.) : x=%f\" x)); 
+  if x > 2. then 2.  else
 
 * AO_coef : float
 

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

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