LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-05-07 07:34:55 +02:00
Code Issues Releases Wiki Activity

minor changes

Browse Source
This commit is contained in:
Emmanuel Giner 2017-03-24 22:31:06 +01:00
parent 2c6cbbc30b
commit 89ae9650aa
9 changed files with 273 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/MRPT/MRPT_Utils.main.irp.f
View File

@ -29,7 +29,7 @@ subroutine routine_3
enddo enddo
enddo enddo
if(save_heff_eigenvectors)then 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 endif
if(N_states.gt.1)then if(N_states.gt.1)then
print*, 'Energy differences : E(i) - E(0)' print*, 'Energy differences : E(i) - E(0)'

7
plugins/MRPT_Utils/EZFIO.cfg
View File

@ -5,3 +5,10 @@ interface: ezfio,provider,ocaml
default: True default: True
[save_heff_eigenvectors]
type: logical
doc: If true, save the eigenvectors of the dressed matrix
interface: ezfio,provider,ocaml
default: False

47
plugins/MRPT_Utils/ezfio_interface.irp.f
View File

@ -1,10 +1,6 @@
! DO NOT MODIFY BY HAND ! DO NOT MODIFY BY HAND
! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py ! 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/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 ] BEGIN_PROVIDER [ logical, do_third_order_1h1p ]
@ -25,12 +21,11 @@ BEGIN_PROVIDER [ logical, do_third_order_1h1p ]
endif endif
END_PROVIDER END_PROVIDER
<<<<<<< HEAD
BEGIN_PROVIDER [ logical, save_heff_eigenvectors ] BEGIN_PROVIDER [ logical, save_heff_eigenvectors ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! If true, you save the eigenvectors of the effective hamiltonian ! If true, save the eigenvectors of the dressed matrix
END_DOC END_DOC
logical :: has logical :: has
@ -45,43 +40,3 @@ BEGIN_PROVIDER [ logical, save_heff_eigenvectors ]
endif endif
END_PROVIDER 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

12
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -181,7 +181,7 @@
accu = 0.d0 accu = 0.d0
do i_state = 1, N_states do i_state = 1, N_states
do i = 1, N_det 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 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) accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
enddo enddo
@ -245,11 +245,11 @@ END_PROVIDER
integer, allocatable :: iorder(:) integer, allocatable :: iorder(:)
! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors ! 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 j=1,min(N_states_diag,N_det)
do i=1,N_det ! do i=1,N_det
CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j) ! CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j)
enddo ! enddo
enddo ! enddo
do j=N_det+1,N_states_diag do j=N_det+1,N_states_diag
do i=1,N_det do i=1,N_det

45
src/AO_Basis/ao_overlap.irp.f
View File

@ -129,3 +129,48 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
END_PROVIDER 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

6
src/Determinants/print_wf.irp.f
View File

@ -40,19 +40,19 @@ subroutine routine
else else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1)) norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif 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 double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,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*,'hmono = ',hmono
print*,'hdouble = ',hdouble print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij print*,'hij = ',hij
else else if (degree == 2)then
print*,'s1',s1 print*,'s1',s1
print*,'h1,p1 = ',h1,p1 print*,'h1,p1 = ',h1,p1
print*,'s2',s2 print*,'s2',s2
print*,'h2,p2 = ',h2,p2 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 endif
print*,'<Ref| H |D_I> = ',hij print*,'<Ref| H |D_I> = ',hij

84
src/MO_Basis/cholesky_mo.irp.f
View File

@ -78,3 +78,87 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align,
enddo enddo
END_PROVIDER 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

122
src/MO_Basis/mos.irp.f
View File

@ -181,24 +181,146 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
allocate ( T(mo_tot_num_align,ao_num) ) allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
! SC
call dgemm('N','N', ao_num, mo_tot_num, ao_num, & call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap,size(ao_overlap,1), & 1.d0, ao_overlap,size(ao_overlap,1), &
mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), &
0.d0, SC, ao_num_align) 0.d0, SC, ao_num_align)
! A.CS
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, & 1.d0, A_mo,LDA_mo, &
SC, size(SC,1), & SC, size(SC,1), &
0.d0, T, mo_tot_num_align) 0.d0, T, mo_tot_num_align)
! SC.A.CS
call dgemm('N','N', ao_num, ao_num, mo_tot_num, & call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC,size(SC,1), & 1.d0, SC,size(SC,1), &
T, mo_tot_num_align, & T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao) 0.d0, A_ao, LDA_ao)
! C(S.A.S)C
! SC.A.CS
deallocate(T,SC) deallocate(T,SC)
end 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) subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none implicit none
BEGIN_DOC BEGIN_DOC

4
src/Utils/LinearAlgebra.irp.f
View File

@ -19,6 +19,10 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
double precision,allocatable :: A_tmp(:,:) double precision,allocatable :: A_tmp(:,:)
allocate (A_tmp(LDA,n)) allocate (A_tmp(LDA,n))
print*, ''
do i = 1, n
print*, A(i,i)
enddo
A_tmp = A A_tmp = A
! Find optimal size for temp arrays ! Find optimal size for temp arrays