diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f index 3b559232..59740877 100644 --- a/plugins/MRPT/MRPT_Utils.main.irp.f +++ b/plugins/MRPT/MRPT_Utils.main.irp.f @@ -29,7 +29,7 @@ subroutine routine_3 enddo enddo if(save_heff_eigenvectors)then - call save_wavefunction_general(N_det_ref,N_states_diag_heff,psi_ref,N_det_ref,CI_dressed_pt2_new_eigenvectors) + call save_wavefunction_general(N_det_ref,N_states_diag,psi_ref,N_det_ref,CI_dressed_pt2_new_eigenvectors) endif if(N_states.gt.1)then print*, 'Energy differences : E(i) - E(0)' diff --git a/plugins/MRPT_Utils/EZFIO.cfg b/plugins/MRPT_Utils/EZFIO.cfg index 2fcc26ad..3342e8bb 100644 --- a/plugins/MRPT_Utils/EZFIO.cfg +++ b/plugins/MRPT_Utils/EZFIO.cfg @@ -5,3 +5,10 @@ interface: ezfio,provider,ocaml default: True +[save_heff_eigenvectors] +type: logical +doc: If true, save the eigenvectors of the dressed matrix +interface: ezfio,provider,ocaml +default: False + + diff --git a/plugins/MRPT_Utils/ezfio_interface.irp.f b/plugins/MRPT_Utils/ezfio_interface.irp.f index 3112b9b6..88abeccc 100644 --- a/plugins/MRPT_Utils/ezfio_interface.irp.f +++ b/plugins/MRPT_Utils/ezfio_interface.irp.f @@ -1,10 +1,6 @@ ! DO NOT MODIFY BY HAND ! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py -<<<<<<< HEAD ! from file /home/giner/qp_fork/quantum_package/src/MRPT_Utils/EZFIO.cfg -======= -! from file /home/scemama/quantum_package/src/MRPT_Utils/EZFIO.cfg ->>>>>>> 4a552cc8fe36ae7c8c86eb714c2f032b44330ea0 BEGIN_PROVIDER [ logical, do_third_order_1h1p ] @@ -25,12 +21,11 @@ BEGIN_PROVIDER [ logical, do_third_order_1h1p ] endif END_PROVIDER -<<<<<<< HEAD BEGIN_PROVIDER [ logical, save_heff_eigenvectors ] implicit none BEGIN_DOC -! If true, you save the eigenvectors of the effective hamiltonian +! If true, save the eigenvectors of the dressed matrix END_DOC logical :: has @@ -45,43 +40,3 @@ BEGIN_PROVIDER [ logical, save_heff_eigenvectors ] endif END_PROVIDER - -BEGIN_PROVIDER [ integer, n_states_diag_heff ] - implicit none - BEGIN_DOC -! Number of eigenvectors obtained with the effective hamiltonian - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrpt_utils_n_states_diag_heff(has) - if (has) then - call ezfio_get_mrpt_utils_n_states_diag_heff(n_states_diag_heff) - else - print *, 'mrpt_utils/n_states_diag_heff not found in EZFIO file' - stop 1 - endif - -END_PROVIDER - -BEGIN_PROVIDER [ logical, pure_state_specific_mrpt2 ] - implicit none - BEGIN_DOC -! If true, diagonalize the dressed matrix for each state and do a state following of the initial states - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrpt_utils_pure_state_specific_mrpt2(has) - if (has) then - call ezfio_get_mrpt_utils_pure_state_specific_mrpt2(pure_state_specific_mrpt2) - else - print *, 'mrpt_utils/pure_state_specific_mrpt2 not found in EZFIO file' - stop 1 - endif - -END_PROVIDER -======= ->>>>>>> 4a552cc8fe36ae7c8c86eb714c2f032b44330ea0 diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index d7b1f0f6..2a424c55 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -181,7 +181,7 @@ accu = 0.d0 do i_state = 1, N_states do i = 1, N_det -! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) + write(*,'(1000(F16.10,x))')delta_ij(i,:,:) do j = i_state, N_det accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) enddo @@ -245,11 +245,11 @@ END_PROVIDER integer, allocatable :: iorder(:) ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors - do j=1,min(N_states_diag,N_det) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j) - enddo - enddo +! do j=1,min(N_states_diag,N_det) +! do i=1,N_det +! CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j) +! enddo +! enddo do j=N_det+1,N_states_diag do i=1,N_det diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index edf48b25..08e57f73 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -129,3 +129,48 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] !$OMP END PARALLEL DO END_PROVIDER +BEGIN_PROVIDER [ double precision, ao_overlap_inv, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Inverse of the overlap matrix + END_DOC + call invert_matrix(ao_overlap, size(ao_overlap,1), ao_num, ao_overlap_inv, size(ao_overlap_inv,1)) +END_PROVIDER + +BEGIN_PROVIDER [double precision, ao_overlap_inv_1_2, (ao_num_align,ao_num)] + implicit none + integer :: i,j,k,l + double precision :: eigvalues(ao_num),eigvectors(ao_num_align, ao_num) + call lapack_diag(eigvalues,eigvectors,ao_overlap,ao_num_align,ao_num) + ao_overlap_inv_1_2 = 0.d0 + double precision :: a_n + do i = 1, ao_num + a_n = 1.d0/dsqrt(eigvalues(i)) + if(a_n.le.1.d-10)cycle + do j = 1, ao_num + do k = 1, ao_num + ao_overlap_inv_1_2(k,j) += eigvectors(k,i) * eigvectors(j,i) * a_n + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, ao_overlap_1_2, (ao_num_align,ao_num)] + implicit none + integer :: i,j,k,l + double precision :: eigvalues(ao_num),eigvectors(ao_num_align, ao_num) + call lapack_diag(eigvalues,eigvectors,ao_overlap,ao_num_align,ao_num) + ao_overlap_1_2 = 0.d0 + double precision :: a_n + do i = 1, ao_num + a_n = dsqrt(eigvalues(i)) + do j = 1, ao_num + do k = 1, ao_num + ao_overlap_1_2(k,j) += eigvectors(k,i) * eigvectors(j,i) * a_n + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/Determinants/print_wf.irp.f b/src/Determinants/print_wf.irp.f index af109e2d..042dc648 100644 --- a/src/Determinants/print_wf.irp.f +++ b/src/Determinants/print_wf.irp.f @@ -40,19 +40,19 @@ subroutine routine else norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1)) endif - print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map) +! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map) double precision :: hmono,hdouble call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble) print*,'hmono = ',hmono print*,'hdouble = ',hdouble print*,'hmono+hdouble = ',hmono+hdouble print*,'hij = ',hij - else + else if (degree == 2)then print*,'s1',s1 print*,'h1,p1 = ',h1,p1 print*,'s2',s2 print*,'h2,p2 = ',h2,p2 - print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) +! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map) endif print*,' = ',hij diff --git a/src/MO_Basis/cholesky_mo.irp.f b/src/MO_Basis/cholesky_mo.irp.f index 97b6abd2..631ffd42 100644 --- a/src/MO_Basis/cholesky_mo.irp.f +++ b/src/MO_Basis/cholesky_mo.irp.f @@ -78,3 +78,87 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, enddo END_PROVIDER +subroutine svd_mo(n,m,P,LDP,C,LDC) + implicit none + BEGIN_DOC +! Singular value decomposition of the AO Density matrix +! +! n : Number of AOs +! +! m : Number of MOs +! +! P(LDP,n) : Density matrix in AO basis +! +! C(LDC,m) : MOs +! + END_DOC + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:), D(:) + + allocate(W(LDC,n),work(2*n),D(n)) + print*, '' + do i = 1, n + print*, P(i,i) + enddo + call svd(P,LDP,C,LDC,D,W,size(W,1),m,n) + double precision :: accu + accu = 0.d0 + print*, 'm',m + do i = 1, m + print*, D(i) + accu += D(i) + enddo + print*,'Sum of D',accu + + deallocate(W,work) +end + +subroutine svd_mo_new(n,m,m_physical,P,LDP,C,LDC) + implicit none + BEGIN_DOC +! Singular value decomposition of the AO Density matrix +! +! n : Number of AOs + +! m : Number of MOs +! +! P(LDP,n) : Density matrix in AO basis +! +! C(LDC,m) : MOs +! +! tol_in : tolerance +! +! rank : Nomber of local MOs (output) +! + END_DOC + integer, intent(in) :: n,m,m_physical, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:), D(:) + + allocate(W(LDC,n),work(2*n),D(n)) + call svd(P,LDP,C,LDC,D,W,size(W,1),m_physical,n) + double precision :: accu + accu = 0.d0 + print*, 'm',m_physical + do i = 1, m_physical + print*, D(i) + accu += D(i) + enddo + print*,'Sum of D',accu + + deallocate(W,work) +end + diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 69abf7b3..1947d856 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -181,24 +181,146 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) allocate ( T(mo_tot_num_align,ao_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T +! SC 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, SC, ao_num_align) +! A.CS call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & 1.d0, A_mo,LDA_mo, & SC, size(SC,1), & 0.d0, T, mo_tot_num_align) +! SC.A.CS call dgemm('N','N', ao_num, ao_num, mo_tot_num, & 1.d0, SC,size(SC,1), & T, mo_tot_num_align, & 0.d0, A_ao, LDA_ao) +! C(S.A.S)C +! SC.A.CS deallocate(T,SC) end + +subroutine mo_to_ao_s_inv_1_2(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis using the S^{-1} matrix + ! S^{-1} C A C^{+} S^{-1} + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo) + double precision, intent(out) :: A_ao(LDA_ao) + double precision, allocatable :: T(:,:), SC_inv_1_2(:,:) + + allocate ( SC_inv_1_2(ao_num_align,mo_tot_num) ) + allocate ( T(mo_tot_num_align,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + +! SC_inv_1_2 = S^{-1}C + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, ao_overlap_inv_1_2,size(ao_overlap_inv_1_2,1), & + mo_coef, size(mo_coef,1), & + 0.d0, SC_inv_1_2, ao_num_align) + +! T = A.(SC_inv_1_2)^{+} + call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & + 1.d0, A_mo,LDA_mo, & + SC_inv_1_2, size(SC_inv_1_2,1), & + 0.d0, T, mo_tot_num_align) + +! SC_inv_1_2.A.CS + call dgemm('N','N', ao_num, ao_num, mo_tot_num, & + 1.d0, SC_inv_1_2,size(SC_inv_1_2,1), & + T, mo_tot_num_align, & + 0.d0, A_ao, LDA_ao) + +! C(S.A.S)C +! SC_inv_1_2.A.CS + deallocate(T,SC_inv_1_2) +end + +subroutine mo_to_ao_s_1_2(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis using the S^{-1} matrix + ! S^{-1} C A C^{+} S^{-1} + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo) + double precision, intent(out) :: A_ao(LDA_ao) + double precision, allocatable :: T(:,:), SC_1_2(:,:) + + allocate ( SC_1_2(ao_num_align,mo_tot_num) ) + allocate ( T(mo_tot_num_align,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + +! SC_1_2 = S^{-1}C + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, ao_overlap_1_2,size(ao_overlap_1_2,1), & + mo_coef, size(mo_coef,1), & + 0.d0, SC_1_2, ao_num_align) + +! T = A.(SC_1_2)^{+} + call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & + 1.d0, A_mo,LDA_mo, & + SC_1_2, size(SC_1_2,1), & + 0.d0, T, mo_tot_num_align) + +! SC_1_2.A.CS + call dgemm('N','N', ao_num, ao_num, mo_tot_num, & + 1.d0, SC_1_2,size(SC_1_2,1), & + T, mo_tot_num_align, & + 0.d0, A_ao, LDA_ao) + +! C(S.A.S)C +! SC_1_2.A.CS + deallocate(T,SC_1_2) +end + + +subroutine mo_to_ao_s_inv(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis using the S^{-1} matrix + ! S^{-1} C A C^{+} S^{-1} + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo) + double precision, intent(out) :: A_ao(LDA_ao) + double precision, allocatable :: T(:,:), SC_inv(:,:) + + allocate ( SC_inv(ao_num_align,mo_tot_num) ) + allocate ( T(mo_tot_num_align,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + +! SC_inv = S^{-1}C + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, ao_overlap_inv,size(ao_overlap_inv,1), & + mo_coef, size(mo_coef,1), & + 0.d0, SC_inv, ao_num_align) + +! T = A.(SC_inv)^{+} + call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & + 1.d0, A_mo,LDA_mo, & + SC_inv, size(SC_inv,1), & + 0.d0, T, mo_tot_num_align) + +! SC_inv.A.CS + call dgemm('N','N', ao_num, ao_num, mo_tot_num, & + 1.d0, SC_inv,size(SC_inv,1), & + T, mo_tot_num_align, & + 0.d0, A_ao, LDA_ao) + +! C(S.A.S)C +! SC_inv.A.CS + deallocate(T,SC_inv) +end + + subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) implicit none BEGIN_DOC diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 44a15ddf..543be36f 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -19,6 +19,10 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) double precision,allocatable :: A_tmp(:,:) allocate (A_tmp(LDA,n)) + print*, '' + do i = 1, n + print*, A(i,i) + enddo A_tmp = A ! Find optimal size for temp arrays