Natural orbitals using SVD

config/ifort.cfg

@@ -31,14 +31,14 @@ OPENMP  : 1          ; Append OpenMP flags
 # -ftz                       : Flushes denormal results to zero
 #
 [OPT]
-FCFLAGS  : -xSSE4.2 -O2 -ip -opt-prefetch -ftz -g
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz -g 
 
 # Profiling flags
 #################
 #
 [PROFILE]
 FC       : -p -g
-FCFLAGS  : -xSSE4.2 -O2 -ip -opt-prefetch -ftz
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz
 
 # Debugging flags
 #################

ocaml/.gitignore

@@ -1,6 +1,7 @@
 _build
 ezfio.ml
 .gitignore
+Git.ml
 Input_auto_generated.ml
 Input_determinants.ml
 Input_hartree_fock.ml

plugins/CAS_SD/cas_s_selected.irp.f

@@ -30,7 +30,7 @@ program full_ci
   endif
 
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
-    call H_apply_CAS_S_selected(pt2, norm_pert, H_pert_diag,  N_st)
+    call H_apply_CAS_S_selected_no_skip(pt2, norm_pert, H_pert_diag,  N_st)
 
     PROVIDE  psi_coef
     PROVIDE  psi_det

scripts/compilation/qp_create_ninja.py

@@ -637,7 +637,7 @@ def ninja_binaries_rule():
     # c m d #
     # ~#~#~ #
 
-    l_cmd = ["cd $module_abs/IRPF90_temp", "ninja $out && [[ -x $out ]] && touch $out"]
+    l_cmd = ["cd $module_abs/IRPF90_temp", "ninja $out && for i in $out ; do [ -x $$i ] && touch $$i ; done"]
 
     # ~#~#~#~#~#~ #
     # s t r i n g #

src/Determinants/density_matrix.irp.f

@@ -10,7 +10,7 @@
    double precision               :: ck, cl, ckl
    double precision               :: phase
    integer                        :: h1,h2,p1,p2,s1,s2, degree
-   integer                        :: exc(0:2,2,2),n_occ_alpha
+   integer                        :: exc(0:2,2,2),n_occ(2)
    double precision, allocatable  :: tmp_a(:,:), tmp_b(:,:)
 
    if(only_single_double_dm)then
@@ -22,7 +22,7 @@
      one_body_dm_mo_beta  = 0.d0
      !$OMP PARALLEL DEFAULT(NONE)                                         &
         !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
-        !$OMP  tmp_a, tmp_b, n_occ_alpha)&
+        !$OMP  tmp_a, tmp_b, n_occ)&
         !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,&
         !$OMP  elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
         !$OMP  mo_tot_num)
@@ -31,8 +31,7 @@
      tmp_b = 0.d0
      !$OMP DO SCHEDULE(dynamic)
      do k=1,N_det
-       call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
+       call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
-       call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
        do m=1,N_states
          ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m)
          do l=1,elec_alpha_num
@@ -182,13 +181,10 @@ subroutine set_natural_mos
  END_DOC
  character*(64) :: label
  double precision, allocatable :: tmp(:,:)
- allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2)))
 
- ! Negation to have the occupied MOs first after the diagonalization
- tmp = one_body_dm_mo
  label = "Natural"
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,-1)
+! call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,label,-1)
- deallocate(tmp)
+ 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)
 
 end
 subroutine save_natural_mos

src/MO_Basis/mo_overlap.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)]
   integer :: i,j,n,l
   double precision :: f
   integer :: lmax
-  lmax = iand(ao_num,-4)
+  lmax = (ao_num/4) * 4
   !$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
   !$OMP  PRIVATE(i,j,n,l) &
   !$OMP  SHARED(mo_overlap,mo_coef,ao_overlap, &

src/MO_Basis/utils.irp.f

@@ -100,6 +100,53 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
   mo_label = label
 end
 
+subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
+  implicit none
+  integer,intent(in)             :: lda,m,n
+  character*(64), intent(in)     :: label
+  double precision, intent(in)   :: matrix(lda,n)
+  
+  integer :: i,j
+  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(output_mo_basis)
+  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_align,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 (output_mo_basis,'(A)'), 'MOs are now **'//trim(label)//'**'
+  write (output_mo_basis,'(A)'), ''
+  write (output_mo_basis,'(A)'), 'Eigenvalues'
+  write (output_mo_basis,'(A)'), '-----------'
+  write (output_mo_basis,'(A)'), ''
+  write (output_mo_basis,'(A)'), '======== ================'
+
+  do i=1,m
+    write (output_mo_basis,'(I8,X,F16.10)'), i,D(i)
+  enddo
+  write (output_mo_basis,'(A)'), '======== ================'
+  write (output_mo_basis,'(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))
+  deallocate(A,mo_coef_new,U,Vt,D)
+  call write_time(output_mo_basis)
+  
+  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

src/Utils/LinearAlgebra.irp.f

@@ -70,9 +70,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   double precision               :: Vt(lda,n)
   double precision               :: D(n)
   double precision               :: S_half(lda,n)
-  double precision,allocatable   :: work(:)
+  !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
-  !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D, work
+  integer                        :: info, i, j, k
-  integer                        :: info, lwork, i, j, k
   
   call svd(overlap,lda,U,ldc,D,Vt,lda,m,n)
 
@@ -82,7 +81,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
 
   !$OMP DO
   do i=1,n
-    if ( D(i) < 1.d-12 ) then
+    if ( D(i) < 1.d-6 ) then
       D(i) = 0.d0
     else
       D(i) = 1.d0/dsqrt(D(i))