Test fock orb

src/BiInts/map_integrals.irp.f

@@ -37,10 +37,10 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
   integer*8, intent(in)          :: i1
   integer*8                      :: i2,i3
   i = 0
-  i2   = ceiling(0.5*(sqrt(8.*real(i1)+1.)-1.))
-  l(1) = ceiling(0.5*(sqrt(8.*real(i2)+1.)-1.))
+  i2   = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0))
+  l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
   i3   = i1 - ishft(i2*i2-i2,-1)
-  k(1) = ceiling(0.5*(sqrt(8.*real(i3)+1.)-1.))
+  k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
   j(1) = i2 - ishft(l(1)*l(1)-l(1),-1)
   i(1) = i3 - ishft(k(1)*k(1)-k(1),-1)
 

src/Dets/save_fock_orb.irp.f

@@ -0,0 +1,45 @@
+program save_fock_orb
+
+  double precision, allocatable :: one_body_dm_ao_alpha(:,:)
+  double precision, allocatable :: one_body_dm_ao_beta (:,:)
+
+  read_wf = .True.
+  touch read_wf
+  print *,  'N_det = ',N_det
+  print *, HF_energy
+  allocate (one_body_dm_ao_alpha(ao_num_align,ao_num),               &
+            one_body_dm_ao_beta (ao_num_align,ao_num))
+
+  call mo_to_ao_no_overlap(one_body_dm_mo_alpha,mo_tot_num_align, &
+                one_body_dm_ao_alpha,ao_num_align)
+  call mo_to_ao_no_overlap(one_body_dm_mo_beta ,mo_tot_num_align, &
+                one_body_dm_ao_beta ,ao_num_align)
+
+  do i=1,mo_tot_num
+   do j=1,mo_tot_num
+     if (abs(fock_matrix_mo(i,j)) > 1.d-10) then
+       print *, i,j, fock_matrix_mo(i,j)
+     endif
+   enddo
+  enddo
+  hf_density_matrix_ao_alpha = one_body_dm_ao_alpha
+  hf_density_matrix_ao_beta  = one_body_dm_ao_beta
+  touch hf_density_matrix_ao_alpha hf_density_matrix_ao_beta
+
+  print *,  '---'
+  do i=1,mo_tot_num
+   do j=1,mo_tot_num
+     if (abs(fock_matrix_mo(i,j)) > 1.d-10) then
+       print *, i,j, fock_matrix_mo(i,j)
+     endif
+   enddo
+  enddo
+  mo_coef = eigenvectors_fock_matrix_mo
+  mo_label = 'CASSCF'
+  TOUCH mo_coef mo_label 
+  print *, HF_energy
+  call save_mos
+
+  deallocate(one_body_dm_ao_alpha,one_body_dm_ao_beta)
+end
+

src/Dets/save_natorb.irp.f

@@ -1,4 +1,6 @@
 program save_natorb
+  read_wf = .True.
+  touch read_wf
   call save_natural_mos
 end
 

src/Hartree_Fock/Fock_matrix.irp.f

@@ -198,11 +198,9 @@ END_PROVIDER
    do i8=0_8,ao_integrals_map%map_size
      n_elements = n_elements_max
      call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
-     if (n_elements == 0) then
-       cycle
-     endif
      do k1=1,n_elements
        call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
+
        do k2=1,8
          if (kk(k2)==0) then
            cycle

src/Hartree_Fock/HF_density_matrix_ao.irp.f

@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Alpha density matrix in the AO basis
+   ! S^-1 x Alpha density matrix in the AO basis x S^-1
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@@ -14,7 +14,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Beta density matrix in the AO basis
+   ! S^-1 Beta density matrix in the AO basis x S^-1
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@@ -27,7 +27,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Density matrix in the AO basis
+   ! S^-1 Density matrix in the AO basis S^-1
    END_DOC
    ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
    if (elec_alpha_num== elec_beta_num) then

src/MOs/mos.irp.f

@@ -86,6 +86,19 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ]
   
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Product S.C where S is the overlap matrix in the AO basis and C the mo_coef matrix.
+ END_DOC
+
+ call dgemm('N','N', ao_num, mo_tot_num, ao_num,                   &
+     1.d0, ao_overlap,size(ao_overlap,1),      &
+     mo_coef, size(mo_coef,1),                                     &
+     0.d0, S_mo_coef, size(S_mo_coef,1))
+
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ]
   implicit none
   BEGIN_DOC

src/Utils/LinearAlgebra.irp.f

@@ -422,3 +422,30 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
   
   deallocate(T,SC)
 end
+
+subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the MO basis to the S^-1 AO basis
+  END_DOC
+  double precision, intent(in)   :: A_mo(LDA_mo)
+  double precision, intent(out)  :: A_ao(LDA_ao)
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  double precision, allocatable  :: T(:,:)
+  
+  allocate ( T(mo_tot_num_align,ao_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
+      1.d0, A_mo,LDA_mo,                                             &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0, T, mo_tot_num_align)
+  
+  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
+      1.d0, mo_coef,size(mo_coef,1),                                 &
+      T, mo_tot_num_align,                                           &
+      0.d0, A_ao, LDA_ao)
+  
+  deallocate(T)
+end
+