diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index 21b91190..bf2026d5 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/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) diff --git a/src/Dets/save_fock_orb.irp.f b/src/Dets/save_fock_orb.irp.f new file mode 100644 index 00000000..67c27db2 --- /dev/null +++ b/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 + diff --git a/src/Dets/save_natorb.irp.f b/src/Dets/save_natorb.irp.f index dd3d7684..e56f9821 100644 --- a/src/Dets/save_natorb.irp.f +++ b/src/Dets/save_natorb.irp.f @@ -1,4 +1,6 @@ program save_natorb + read_wf = .True. + touch read_wf call save_natural_mos end diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index c2c98630..e722d8bf 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/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 diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index e85e4a6c..e8585f59 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/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 diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index 1a0cea77..f1158530 100644 --- a/src/MOs/mos.irp.f +++ b/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 diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 8759bf0f..9f2336db 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/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 +