LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
Code Issues Releases Wiki Activity

two bod is ok

Browse Source
This commit is contained in:
Emmanuel Giner 2016-04-17 22:25:25 +02:00
parent 80cf1472ca
commit 74f465be90
2 changed files with 141 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
src/Determinants/two_body_dm_map.irp.f
View File

@ -8,6 +8,7 @@ BEGIN_PROVIDER [ type(map_type), two_body_dm_ab_map ]
END_DOC END_DOC
integer(key_kind) :: key_max integer(key_kind) :: key_max
integer(map_size_kind) :: sze integer(map_size_kind) :: sze
use map_module
call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max) call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
sze = key_max sze = key_max
call map_init(two_body_dm_ab_map,sze) call map_init(two_body_dm_ab_map,sze)
@ -129,51 +130,56 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree>2)cycle if(degree>2)cycle
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
contrib = psi_coef(i,1) * psi_coef(j,1) * phase
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************
if(s1==s2)cycle ! Only the alpha/beta two body density matrix if(s1==s2)cycle ! Only the alpha/beta two body density matrix
! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J ! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
if(h1>p1)cycle
if(h2>p2)cycle
! if(s1.ne.1)cycle
n_elements += 1 n_elements += 1
contrib = psi_coef(i,1) * psi_coef(j,1) * phase
buffer_value(n_elements) = contrib buffer_value(n_elements) = contrib
!DEC$ FORCEINLINE !DEC$ FORCEINLINE
! call mo_bielec_integrals_index(h1,p1,h2,p2,buffer_i(n_elements))
call mo_bielec_integrals_index(h1,h2,p1,p2,buffer_i(n_elements)) call mo_bielec_integrals_index(h1,h2,p1,p2,buffer_i(n_elements))
if (n_elements == size_buffer) then ! if (n_elements == size_buffer) then
call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& ! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind)) ! real(mo_integrals_threshold,integral_kind))
n_elements = 0 ! n_elements = 0
endif ! endif
else ! case of the SINGLE EXCITATIONS *************************************************** else ! case of the SINGLE EXCITATIONS ***************************************************
cycle
if(s1==1)then ! Mono alpha : ! if(s1==1)then ! Mono alpha :
do k = 1, elec_beta_num ! do k = 1, elec_beta_num
m = occ(k,2) ! m = occ(k,2)
n_elements += 1 ! n_elements += 1
buffer_value(n_elements) = contrib ! buffer_value(n_elements) = contrib
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J ! ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements)) ! call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
if (n_elements == size_buffer) then ! if (n_elements == size_buffer) then
call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& ! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind)) ! real(mo_integrals_threshold,integral_kind))
n_elements = 0 ! n_elements = 0
endif ! endif
enddo ! enddo
else ! Mono Beta : ! else ! Mono Beta :
do k = 1, elec_alpha_num ! do k = 1, elec_alpha_num
m = occ(k,1) ! m = occ(k,1)
n_elements += 1 ! n_elements += 1
buffer_value(n_elements) = contrib ! buffer_value(n_elements) = contrib
! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J ! ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements)) ! call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
if (n_elements == size_buffer) then ! if (n_elements == size_buffer) then
call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& ! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind)) ! real(mo_integrals_threshold,integral_kind))
n_elements = 0 ! n_elements = 0
endif ! endif
enddo ! enddo
endif ! endif
endif endif
enddo enddo
@ -181,6 +187,9 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
print*,'n_elements = ',n_elements print*,'n_elements = ',n_elements
call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind)) real(mo_integrals_threshold,integral_kind))
call map_unique(two_body_dm_ab_map)
deallocate(buffer_i,buffer_value)
end end
@ -192,8 +201,8 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)]
double precision :: contrib double precision :: contrib
BEGIN_DOC BEGIN_DOC
! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi> ! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi>
END_DOC END_DOC
two_body_dm_ab_diag = 0.d0 two_body_dm_ab_diag = 0.d0
do i = 1, N_det ! i == |I> do i = 1, N_det ! i == |I>
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
@ -202,9 +211,83 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)]
k = occ(j,2) k = occ(j,2)
do l = 1, elec_beta_num do l = 1, elec_beta_num
m = occ(l,1) m = occ(l,1)
two_body_dm_ab_diag(k,m) += contrib two_body_dm_ab_diag(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag(m,k) += 0.5d0 * contrib
enddo enddo
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
integer :: i,j,k,l,m
integer :: degree
PROVIDE mo_coef psi_coef psi_det
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
double precision :: phase
double precision :: contrib
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
two_body_dm_ab_big_array = 0.d0
BEGIN_DOC
! The alpha-beta energy can be computed thanks to
! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array(h1,p1,h2,p2) * (h1p1|h2p2)
END_DOC
do i = 1, N_det ! i == |I>
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do j = i+1, N_det ! j == <J|
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree>2)cycle
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase
if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************
if(s1==s2)cycle ! Only the alpha/beta two body density matrix
! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
else if(degree==1)then! case of the SINGLE EXCITATIONS ***************************************************
if(s1==1)then ! Mono alpha :
do k = 1, elec_beta_num
m = occ(k,2)
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
else ! Mono Beta :
do k = 1, elec_alpha_num
m = occ(k,1)
! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
endif
endif
enddo
enddo
print*,'Big array for density matrix provided !'
END_PROVIDER
subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contrib,h1,p1,h2,p2)
implicit none
integer, intent(in) :: h1,p1,h2,p2
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4)
double precision :: contrib
big_array(h1,p1,h2,p2) += 1.d0 * contrib
big_array(p1,h1,h2,p2) += 1.d0 * contrib
big_array(h1,p1,p2,h2) += 1.d0 * contrib
big_array(p1,h1,p2,h2) += 1.d0 * contrib
!big_array(h2,p2,h1,p1) += 1.d0 * contrib
!big_array(p2,h2,h1,p1) += 1.d0 * contrib
!if(p2.ne.h2)then
!big_array(h2,p2,p1,h1) += 1.d0 * contrib
!big_array(p2,h2,p1,h1) += 1.d0 * contrib
!endif
end

23
src/MO_Basis/utils.irp.f
View File

@ -268,3 +268,26 @@ subroutine mo_sort_by_observable(observable,label)
end end
subroutine give_all_mos_at_r(r,mos_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
call give_specific_mos_at_r(r,mos_array, mo_coef)
end
subroutine give_specific_mos_at_r(r,mos_array, mo_coef_specific)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(in) :: mo_coef_specific(ao_num_align, mo_tot_num)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision :: aos_array(ao_num),accu
integer :: i,j
call give_all_aos_at_r(r,aos_array)
do i = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
accu += mo_coef_specific(j,i) * aos_array(j)
enddo
mos_array(i) = accu
enddo
end