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

merging with new way

Browse Source
This commit is contained in:
Emmanuel Giner 2016-09-10 17:00:34 +02:00
commit 8ad7a5c82f
8 changed files with 339 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
# 0 : Deactivate # 0 : Deactivate
# #
[OPTION] [OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py CACHE : 1 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags OPENMP : 1 ; Append OpenMP flags

3
plugins/MRPT_Utils/excitations_cas.irp.f
View File

@ -370,9 +370,6 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe
do j = 1, ndet do j = 1, ndet
if(psi_coef_tmp(j)==0.d0)cycle if(psi_coef_tmp(j)==0.d0)cycle
call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
! call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij_bis)
! print*, 'i,j',i,j
! print*, hij
accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij
enddo enddo
enddo enddo

10
plugins/MRPT_Utils/fock_like_operators.irp.f
View File

@ -11,8 +11,6 @@
accu = 0.d0 accu = 0.d0
i_inact_core_orb = list_core_inact(i) i_inact_core_orb = list_core_inact(i)
do j = 1, n_core_inact_orb do j = 1, n_core_inact_orb
! do j = 1, elec_alpha_num
! j_inact_core_orb = j
j_inact_core_orb = list_core_inact(j) j_inact_core_orb = list_core_inact(j)
accu += 2.d0 * mo_bielec_integral_jj(i_inact_core_orb,j_inact_core_orb) & accu += 2.d0 * mo_bielec_integral_jj(i_inact_core_orb,j_inact_core_orb) &
- mo_bielec_integral_jj_exchange(i_inact_core_orb,j_inact_core_orb) - mo_bielec_integral_jj_exchange(i_inact_core_orb,j_inact_core_orb)
@ -84,8 +82,8 @@
accu_exchange(2) += 2.d0 * nb * exchange accu_exchange(2) += 2.d0 * nb * exchange
enddo enddo
enddo enddo
fock_core_inactive_from_act(i_inact_core_orb,1,i_state) = accu_coulomb + accu_exchange(1) fock_core_inactive_from_act(i_inact_core_orb,1,i_state) = accu_coulomb - accu_exchange(1)
fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb + accu_exchange(2) fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb - accu_exchange(2)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
@ -131,8 +129,8 @@
accu_exchange(2) += 2.d0 * nb * exchange accu_exchange(2) += 2.d0 * nb * exchange
enddo enddo
enddo enddo
fock_virt_from_act(i_virt_orb,1,i_state) = accu_coulomb + accu_exchange(1) fock_virt_from_act(i_virt_orb,1,i_state) = accu_coulomb - accu_exchange(1)
fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb + accu_exchange(2) fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb - accu_exchange(2)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER

15
plugins/MRPT_Utils/mrpt_dress.irp.f
View File

@ -44,6 +44,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
integer :: N_miniList, leng integer :: N_miniList, leng
double precision :: delta_e(N_states),hij_tmp double precision :: delta_e(N_states),hij_tmp
integer :: index_i,index_j integer :: index_i,index_j
double precision :: phase_array(N_det),phase
integer :: exc(0:2,2,2),degree
leng = max(N_det_generators, N_det) leng = max(N_det_generators, N_det)
@ -74,11 +76,14 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
! double precision :: ihpsi0,coef_pert ! double precision :: ihpsi0,coef_pert
! ihpsi0 = 0.d0 ! ihpsi0 = 0.d0
! coef_pert = 0.d0 ! coef_pert = 0.d0
phase_array =0.d0
do i = 1,idx_alpha(0) do i = 1,idx_alpha(0)
index_i = idx_alpha(i) index_i = idx_alpha(i)
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha)
hij_array(index_i) = hialpha hij_array(index_i) = hialpha
call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
! phase_array(index_i) = phase
do i_state = 1,N_states do i_state = 1,N_states
delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
enddo enddo
@ -90,6 +95,16 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
call omp_set_lock( psi_ref_bis_lock(index_i) ) call omp_set_lock( psi_ref_bis_lock(index_i) )
do j = 1, idx_alpha(0) do j = 1, idx_alpha(0)
index_j = idx_alpha(j) index_j = idx_alpha(j)
! call get_excitation(psi_det(1,1,index_i),psi_det(1,1,index_i),exc,degree,phase,N_int)
! if(index_j.ne.index_i)then
! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then
! print*, phase_array(index_j) , phase_array(index_i) ,phase
! call debug_det(psi_det(1,1,index_i),N_int)
! call debug_det(psi_det(1,1,index_j),N_int)
! call debug_det(tq(1,1,i_alpha),N_int)
! stop
! endif
! endif
do i_state=1,N_states do i_state=1,N_states
! standard dressing first order ! standard dressing first order
delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state) delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state)

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

@ -130,19 +130,30 @@
print*, '2h1p = ',accu print*, '2h1p = ',accu
! 2h2p ! 2h2p
delta_ij_tmp = 0.d0 !delta_ij_tmp = 0.d0
call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) !call H_apply_mrpt_2h2p(delta_ij_tmp,N_det)
accu = 0.d0 !accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_2h2p(i_state) = accu(i_state)
!enddo
!print*, '2h2p = ',accu
double precision :: contrib_2h2p(N_states)
call give_2h2p(contrib_2h2p)
do i_state = 1, N_states do i_state = 1, N_states
do i = 1, N_det do i = 1, N_det
do j = 1, N_det delta_ij(i,i,i_state) += contrib_2h2p(i_state)
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo enddo
second_order_pt_new_2h2p(i_state) = accu(i_state) second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state)
enddo enddo
print*, '2h2p = ',accu print*, '2h2p = ',contrib_2h2p(1)
! total ! total
accu = 0.d0 accu = 0.d0

192
plugins/MRPT_Utils/new_way.irp.f Normal file
View File

@ -0,0 +1,192 @@
subroutine give_2h1p_contrib(matrix_2h1p)
use bitmasks
implicit none
double precision , intent(inout) :: matrix_2h1p(N_det,N_det,*)
integer :: i,j,r,a,b
integer :: iorb, jorb, rorb, aorb, borb
integer :: ispin,jspin
integer :: idet,jdet
integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2)
double precision :: perturb_dets_phase(n_act_orb,2,2)
double precision :: perturb_dets_hij(n_act_orb,2,2)
double precision :: coef_perturb_from_idet(n_act_orb,2,2,N_states)
integer :: inint
integer :: elec_num_tab_local(2),acu_elec
integer(bit_kind) :: det_tmp(N_int,2)
integer :: exc(0:2,2,2)
integer :: accu_elec
double precision :: get_mo_bielec_integral_schwartz
double precision :: active_int(n_act_orb,2)
double precision :: hij,phase
!matrix_2h1p = 0.d0
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
do j = 1, n_inact_orb ! Second inactive
jorb = list_inact(j)
do r = 1, n_virt_orb ! First virtual
rorb = list_virt(r)
! take all the integral you will need for i,j,r fixed
do a = 1, n_act_orb
aorb = list_act(a)
active_int(a,1) = get_mo_bielec_integral_schwartz(iorb,jorb,rorb,aorb,mo_integrals_map) ! direct
active_int(a,2) = get_mo_bielec_integral_schwartz(iorb,jorb,aorb,rorb,mo_integrals_map) ! exchange
enddo
integer :: degree(N_det)
integer :: idx(0:N_det)
double precision :: delta_e(n_act_orb,2,N_states)
integer :: istate
integer :: index_orb_act_mono(N_det,3)
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a)
if(ispin == jspin .and. iorb.le.jorb)cycle ! condition not to double count
do a = 1, n_act_orb ! First active
aorb = list_act(a)
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
! Do the excitation inactive -- > active
call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin
call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
enddo
if(accu_elec .ne. elec_num_tab_local(jspin))then
perturb_dets_phase(a,jspin,ispin) = 0.0
perturb_dets_hij(a,jspin,ispin) = 0.d0
do istate = 1, N_states
coef_perturb_from_idet(a,jspin,ispin,istate) = 0.d0
enddo
cycle
endif
do inint = 1, N_int
perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_creat(a,jspin,istate) &
- fock_virt_total_spin_trace(rorb,istate) &
+ fock_core_inactive_total_spin_trace(iorb,istate) &
+ fock_core_inactive_total_spin_trace(jorb,istate)
enddo
if(ispin == jspin)then
perturb_dets_hij(a,jspin,ispin) = phase * (active_int(a,2) - active_int(a,1) )
else
perturb_dets_hij(a,jspin,ispin) = phase * active_int(a,1)
endif
!!!!!!!!!!!!!!!!!!!!!1 Computation of the coefficient at first order coming from idet
!!!!!!!!!!!!!!!!!!!!! for the excitation (i,j)(ispin,jspin) ---> (r,a)(ispin,jspin)
do istate = 1, N_states
coef_perturb_from_idet(a,jspin,ispin,istate) = perturb_dets_hij(a,jspin,ispin) / delta_e(a,jspin,istate)
enddo
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a_{b} a^{\dagger}_a | Idet>
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,1)) !!! a^{\dagger}_a
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a_{b}
index_orb_act_mono(idx(jdet),3) = 1
else
! Mono beta
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,2)) !!! a_{b}
index_orb_act_mono(idx(jdet),3) = 2
endif
else
index_orb_act_mono(idx(jdet),1) = -1
endif
enddo
integer :: kspin
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator
! are connected by the presence of the perturbers determinants |det_tmp>
aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb}
borb = index_orb_act_mono(idx(jdet),2) ! a_{borb}
kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation
! the determinants Idet and Jdet interact throw the following operator
! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet >
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin)
enddo
double precision :: hja
! you determine the interaction between the excited determinant and the other parent | Jdet >
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet >
! hja = < det_tmp | H | Jdet >
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,2) - active_int(borb,1) )
else
hja = phase * active_int(borb,1)
endif
do istate = 1, N_states
matrix_2h1p(idx(jdet),idet,istate) += hja * coef_perturb_from_idet(aorb,kspin,ispin,istate)
enddo
enddo ! ispin
else
! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations
!
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
do ispin = 1, 2
do kspin = 1, 2
if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
do a = 1, n_act_orb ! First active
do istate = 1, N_states
matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate) * perturb_dets_hij(a,kspin,ispin)
enddo
enddo
enddo
enddo
endif
enddo
enddo
enddo
enddo
enddo
end

6
plugins/MRPT_Utils/psi_active_prov.irp.f
View File

@ -183,6 +183,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
double precision :: delta_e_inactive(N_states) double precision :: delta_e_inactive(N_states)
integer :: i_hole_inact integer :: i_hole_inact
call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list)
delta_e_inactive = 0.d0 delta_e_inactive = 0.d0
do i = 1, n_holes_spin(1) do i = 1, n_holes_spin(1)
@ -432,7 +434,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
do i_state = 1, n_states do i_state = 1, n_states
delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state) delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state)
enddo enddo
end end

103
src/Determinants/slater_rules.irp.f
View File

@ -749,6 +749,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
exc(1,2,2) ,mo_integrals_map) exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then else if (exc(0,1,1) == 2) then
! Double alpha ! Double alpha
print*,'phase hij = ',phase
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral_schwartz( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
@ -759,6 +760,17 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
exc(2,1,1), & exc(2,1,1), &
exc(2,2,1), & exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) ) exc(1,2,1) ,mo_integrals_map) )
print*,get_mo_bielec_integral_schwartz( &
exc(1,1,1), &
exc(2,1,1), &
exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map)
print*,get_mo_bielec_integral_schwartz( &
exc(1,1,1), &
exc(2,1,1), &
exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map)
else if (exc(0,1,2) == 2) then else if (exc(0,1,2) == 2) then
! Double beta ! Double beta
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral_schwartz( &
@ -1225,6 +1237,97 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a
print*,'------' print*,'------'
end end
subroutine get_excitation_degree_vector_mono(key1,key2,degree,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Applies get_excitation_degree to an array of determinants and return only the mono excitations
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: degree(sze)
integer, intent(out) :: idx(0:sze)
integer :: i,l,d,m
ASSERT (Nint > 0)
ASSERT (sze > 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (d > 2) then
cycle
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (d > 2) then
cycle
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (d > 2) then
cycle
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = 0
!DIR$ LOOP COUNT MIN(4)
do m=1,Nint
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
+ popcnt(xor( key1(m,2,i), key2(m,2)))
enddo
if (d > 2) then
cycle
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)