mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
experimental mrcepa0/mrsc2
This commit is contained in:
Yann Garniron
parent
f65aae9538
commit
7f583946c2
1
ocaml/.gitignore
vendored
1
ocaml/.gitignore
vendored
@ -4,7 +4,6 @@ ezfio.ml
|
||||
Git.ml
|
||||
Input_auto_generated.ml
|
||||
Input_determinants.ml
|
||||
Input_foboci.ml
|
||||
Input_hartree_fock.ml
|
||||
Input_integrals_bielec.ml
|
||||
Input_perturbation.ml
|
||||
|
||||
@ -17,13 +17,12 @@ type keyword =
|
||||
| Electrons
|
||||
| Mo_basis
|
||||
| Nuclei
|
||||
| Hartree_fock
|
||||
| Pseudo
|
||||
| Integrals_bielec
|
||||
| Perturbation
|
||||
| Properties
|
||||
| Foboci
|
||||
| Determinants
|
||||
| Integrals_bielec
|
||||
| Pseudo
|
||||
| Perturbation
|
||||
| Hartree_fock
|
||||
| Properties
|
||||
;;
|
||||
|
||||
|
||||
@ -33,13 +32,12 @@ let keyword_to_string = function
|
||||
| Electrons -> "Electrons"
|
||||
| Mo_basis -> "MO basis"
|
||||
| Nuclei -> "Molecule"
|
||||
| Hartree_fock -> "Hartree_fock"
|
||||
| Pseudo -> "Pseudo"
|
||||
| Integrals_bielec -> "Integrals_bielec"
|
||||
| Perturbation -> "Perturbation"
|
||||
| Properties -> "Properties"
|
||||
| Foboci -> "Foboci"
|
||||
| Determinants -> "Determinants"
|
||||
| Integrals_bielec -> "Integrals_bielec"
|
||||
| Pseudo -> "Pseudo"
|
||||
| Perturbation -> "Perturbation"
|
||||
| Hartree_fock -> "Hartree_fock"
|
||||
| Properties -> "Properties"
|
||||
;;
|
||||
|
||||
|
||||
@ -88,20 +86,18 @@ let get s =
|
||||
f Ao_basis.(read, to_rst)
|
||||
| Determinants_by_hand ->
|
||||
f Determinants_by_hand.(read, to_rst)
|
||||
| Hartree_fock ->
|
||||
f Hartree_fock.(read, to_rst)
|
||||
| Pseudo ->
|
||||
f Pseudo.(read, to_rst)
|
||||
| Integrals_bielec ->
|
||||
f Integrals_bielec.(read, to_rst)
|
||||
| Perturbation ->
|
||||
f Perturbation.(read, to_rst)
|
||||
| Properties ->
|
||||
f Properties.(read, to_rst)
|
||||
| Foboci ->
|
||||
f Foboci.(read, to_rst)
|
||||
| Determinants ->
|
||||
f Determinants.(read, to_rst)
|
||||
| Integrals_bielec ->
|
||||
f Integrals_bielec.(read, to_rst)
|
||||
| Pseudo ->
|
||||
f Pseudo.(read, to_rst)
|
||||
| Perturbation ->
|
||||
f Perturbation.(read, to_rst)
|
||||
| Hartree_fock ->
|
||||
f Hartree_fock.(read, to_rst)
|
||||
| Properties ->
|
||||
f Properties.(read, to_rst)
|
||||
end
|
||||
with
|
||||
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
|
||||
@ -139,13 +135,12 @@ let set str s =
|
||||
in
|
||||
let open Input in
|
||||
match s with
|
||||
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
|
||||
| Pseudo -> write Pseudo.(of_rst, write) s
|
||||
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
|
||||
| Perturbation -> write Perturbation.(of_rst, write) s
|
||||
| Properties -> write Properties.(of_rst, write) s
|
||||
| Foboci -> write Foboci.(of_rst, write) s
|
||||
| Determinants -> write Determinants.(of_rst, write) s
|
||||
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
|
||||
| Pseudo -> write Pseudo.(of_rst, write) s
|
||||
| Perturbation -> write Perturbation.(of_rst, write) s
|
||||
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
|
||||
| Properties -> write Properties.(of_rst, write) s
|
||||
| Electrons -> write Electrons.(of_rst, write) s
|
||||
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
|
||||
| Nuclei -> write Nuclei.(of_rst, write) s
|
||||
@ -193,13 +188,12 @@ let run check_only ezfio_filename =
|
||||
Nuclei ;
|
||||
Ao_basis;
|
||||
Electrons ;
|
||||
Hartree_fock ;
|
||||
Pseudo ;
|
||||
Integrals_bielec ;
|
||||
Perturbation ;
|
||||
Properties ;
|
||||
Foboci ;
|
||||
Determinants ;
|
||||
Integrals_bielec ;
|
||||
Pseudo ;
|
||||
Perturbation ;
|
||||
Hartree_fock ;
|
||||
Properties ;
|
||||
Mo_basis;
|
||||
Determinants_by_hand ;
|
||||
]
|
||||
|
||||
@ -359,7 +359,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
|
||||
y, &
|
||||
lambda &
|
||||
)
|
||||
abort_here = abort_all
|
||||
end
|
||||
|
||||
|
||||
|
||||
@ -1,3 +1,4 @@
|
||||
|
||||
BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ]
|
||||
END_PROVIDER
|
||||
|
||||
@ -47,6 +48,7 @@
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
i_pert = 0
|
||||
if( i_pert == 1)then
|
||||
pert_determinants(k,i) = i_pert
|
||||
endif
|
||||
@ -110,7 +112,6 @@ END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref,N_det_ref,N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Dressing matrix in N_det basis
|
||||
@ -118,7 +119,6 @@ END_PROVIDER
|
||||
integer :: i,j,m
|
||||
delta_ij = 0.d0
|
||||
delta_ii = 0.d0
|
||||
|
||||
call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref)
|
||||
double precision :: max_delta
|
||||
double precision :: accu
|
||||
@ -135,6 +135,7 @@ END_PROVIDER
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
!stop "movais delta"
|
||||
print*,''
|
||||
print*,''
|
||||
print*,'<psi| Delta H |psi> = ',accu
|
||||
@ -159,17 +160,6 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
|
||||
h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!!!!!!!!!!
|
||||
do ii = 1, N_det_ref
|
||||
do jj = 1, N_det_ref
|
||||
i = idx_ref(ii)
|
||||
j = idx_ref(jj)
|
||||
h_matrix_dressed(i,j,istate) += delta_cas(ii,jj,istate)
|
||||
h_matrix_dressed(j,i,istate) += delta_cas(ii,jj,istate)
|
||||
end do
|
||||
end do
|
||||
!!!!!!!!!!!!!
|
||||
do ii = 1, N_det_ref
|
||||
i =idx_ref(ii)
|
||||
h_matrix_dressed(i,i,istate) += delta_ii(ii,istate)
|
||||
@ -272,11 +262,14 @@ subroutine diagonalize_CI_dressed
|
||||
! eigenstates of the CI matrix
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
double precision, parameter :: speed = 1d0
|
||||
|
||||
do j=1,N_states_diag
|
||||
do i=1,N_det
|
||||
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
|
||||
psi_coef(i,j) = CI_eigenvectors_dressed(i,j) * speed + psi_coef(i,j) * (1d0 - speed)
|
||||
enddo
|
||||
enddo
|
||||
SOFT_TOUCH psi_coef
|
||||
|
||||
end
|
||||
|
||||
|
||||
@ -1,25 +1,58 @@
|
||||
use bitmasks
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
! delta_ij(:,:,:) = delta_ij_old(:,:,:)
|
||||
! delta_ii(:,:) = delta_ii_old(:,:)
|
||||
delta_ij(:,:,:) = delta_mrcepa0_ij(:,:,:)! - delta_sub_ij(:,:,:)
|
||||
delta_ii(:,:)= delta_mrcepa0_ii(:,:)! - delta_sub_ii(:,:)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ]
|
||||
&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ]
|
||||
&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_det_non_ref) ]
|
||||
&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_det_ref) ]
|
||||
&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (2) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), active_sorb(2)
|
||||
integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), nonactive_sorb(2)
|
||||
integer i, II, j, k
|
||||
logical, external :: detEq
|
||||
|
||||
print *, "provide cepa0"
|
||||
active_sorb = (/b'001100000000', b'001100000000'/)
|
||||
active_sorb(:) = 0_8
|
||||
nonactive_sorb(:) = not(0_8)
|
||||
|
||||
if(N_det_ref > 1) then
|
||||
do i=1, N_det_ref
|
||||
active_sorb(1) = ior(psi_ref(1,1,i), active_sorb(1))
|
||||
active_sorb(2) = ior(psi_ref(1,2,i), active_sorb(2))
|
||||
nonactive_sorb(1) = iand(psi_ref(1,1,i), nonactive_sorb(1))
|
||||
nonactive_sorb(2) = iand(psi_ref(1,2,i), nonactive_sorb(2))
|
||||
end do
|
||||
active_sorb(1) = iand(active_sorb(1), not(nonactive_sorb(1)))
|
||||
active_sorb(2) = iand(active_sorb(2), not(nonactive_sorb(2)))
|
||||
end if
|
||||
|
||||
do i=1, N_det_non_ref
|
||||
det_noactive(1,1,i) = iand(psi_non_ref(1,1,i), not(active_sorb(1)))
|
||||
det_noactive(1,2,i) = iand(psi_non_ref(1,2,i), not(active_sorb(2)))
|
||||
end do
|
||||
|
||||
call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int)
|
||||
|
||||
! do i=1, N_det_non_ref
|
||||
! print "(B30,B30)", det_noactive(1,1,i), det_noactive(1,2,i)
|
||||
! end do
|
||||
! stop
|
||||
do i=1,N_det_ref
|
||||
det_ref_active(i) = iand(psi_ref(1,1,i), active_sorb(1))
|
||||
det_ref_active(i) = det_ref_active(i) + iand(psi_ref(1,2,i), active_sorb(2)) * 2_8**32_8
|
||||
@ -39,8 +72,8 @@ use bitmasks
|
||||
cepa0_shortcut(cepa0_shortcut(0)) = i
|
||||
end if
|
||||
end do
|
||||
cepa0_shortcut(0) += 1
|
||||
cepa0_shortcut(cepa0_shortcut(0)) = N_det_non_ref+1
|
||||
!cepa0_shortcut(0) += 1
|
||||
cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -50,24 +83,22 @@ END_PROVIDER
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
double precision :: Hjk, Hki
|
||||
double precision :: Hjk, Hki, Hij, mat(2,2)
|
||||
integer i_state
|
||||
|
||||
provide lambda_mrcc
|
||||
i_state = 1
|
||||
|
||||
do i=1,N_det_ref
|
||||
do j=1,i
|
||||
delta_cas(i,j,i_state) = 0d0
|
||||
do k=1,N_det_non_ref
|
||||
call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
|
||||
call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,k),N_int,Hki)
|
||||
call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
|
||||
delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k)
|
||||
end do
|
||||
delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
|
||||
end do
|
||||
end do
|
||||
|
||||
print *, "mrcepa0_cas_dressing", delta_cas(:,:,1)
|
||||
END_PROVIDER
|
||||
|
||||
logical function detEq(a,b,Nint)
|
||||
@ -88,8 +119,9 @@ end function
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
@ -99,14 +131,14 @@ end function
|
||||
double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
|
||||
double precision :: contrib
|
||||
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
|
||||
integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle
|
||||
integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive
|
||||
integer, allocatable :: idx_sorted_bit(:)
|
||||
integer, external :: get_index_in_psi_det_sorted_bit
|
||||
logical, external :: is_in_wavefunction
|
||||
|
||||
integer :: II, blok
|
||||
|
||||
provide det_cepa0_active
|
||||
provide det_cepa0_active delta_cas lambda_mrcc
|
||||
|
||||
if(N_int /= 1) then
|
||||
print *, "mrcepa0 experimental N_int==1"
|
||||
@ -114,11 +146,11 @@ end function
|
||||
end if
|
||||
|
||||
i_state = 1
|
||||
delta_ii(:,:) = 0
|
||||
delta_ij(:,:,:) = 0
|
||||
delta_mrcepa0_ii(:,:) = 0d0
|
||||
delta_mrcepa0_ij(:,:,:) = 0d0
|
||||
|
||||
! do i=1,N_det_ref
|
||||
! delta_ii(i,i_state) = delta_cas(i,i)
|
||||
! delta_ii(i,i_state) = delta_cas(i,i,i_state)
|
||||
! end do
|
||||
|
||||
provide mo_bielec_integrals_in_map
|
||||
@ -128,43 +160,193 @@ end function
|
||||
do i=1,N_det_non_ref
|
||||
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
||||
enddo
|
||||
|
||||
!sd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii)
|
||||
!- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii)
|
||||
do blok=1,cepa0_shortcut(0)
|
||||
do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
||||
do II=1,N_det_ref
|
||||
|
||||
made_hole = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
|
||||
made_particle = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
|
||||
|
||||
if(popcnt(made_hole) + popcnt(made_particle) > 2) cycle
|
||||
made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
|
||||
made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
|
||||
call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int)
|
||||
if (degree > 2 .or. popcnt(made_hole) * popcnt(made_particle) /= degree*2) cycle
|
||||
|
||||
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
||||
if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP"
|
||||
!do k=1,N_det_non_ref
|
||||
if(iand(made_hole, det_cepa0_active(k)) /= 0) cycle
|
||||
if(iand(made_particle, det_cepa0_active(k)) /= made_particle) cycle
|
||||
myActive = xor(det_cepa0_active(k), made_hole)
|
||||
myActive = xor(myActive, made_particle)
|
||||
if(i==k .and. myActive /= det_ref_active(II)) stop "AAAA"
|
||||
!if(i==k) print *, "i=k"
|
||||
do J=1,N_det_ref
|
||||
if(iand(made_hole, det_ref_active(J)) /= made_hole) cycle
|
||||
if(iand(made_particle, det_ref_active(J)) /= 0) cycle
|
||||
call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,det_cepa0_idx(k)),N_int,Hki)
|
||||
if(det_ref_active(J) /= myActive) cycle
|
||||
|
||||
!!!!!
|
||||
call get_excitation_degree(psi_ref(1,1,J),psi_non_ref(1,1,det_cepa0_idx(k)),degree,N_int)
|
||||
if(degree > 2) stop "BBBB"
|
||||
!!!!!!!!!
|
||||
! if(i/=k .and. popcnt(made_hole) /= popcnt(made_particle)) then
|
||||
! print *, "=================", made_hole, made_particle
|
||||
! call debug_det(psi_ref(1,1,II),N_int)
|
||||
! call debug_det(psi_non_ref(1,1,det_cepa0_idx(i)),N_int)
|
||||
! call debug_det(psi_ref(1,1,J),N_int)
|
||||
! call debug_det(psi_non_ref(1,1,det_cepa0_idx(k)),N_int)
|
||||
! print *, "================="
|
||||
! end if
|
||||
|
||||
call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,II),N_int,Hki)
|
||||
|
||||
contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)
|
||||
delta_ij(II, det_cepa0_idx(i), i_state) += contrib
|
||||
delta_mrcepa0_ij(II, det_cepa0_idx(i), i_state) += contrib
|
||||
!
|
||||
if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
|
||||
!-$OMP CRITICAL
|
||||
delta_mrcepa0_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
|
||||
!-$OMP END CRITICAL
|
||||
endif
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!- qs $OMP END PARALLEL DO
|
||||
!print *, "MMMMMMMMMM ", delta_cas(2,2,i_state) , delta_ii(2,i_state)
|
||||
! do i=1,N_det_ref
|
||||
! delta_cas(i,i,i_state) += delta_ii(i,i_state)
|
||||
! end do
|
||||
|
||||
! if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
|
||||
! !qs$OMP CRITICAL
|
||||
! delta_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
|
||||
! !dsd$OMP END CRITICAL
|
||||
! endif
|
||||
deallocate(idx_sorted_bit)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni
|
||||
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
|
||||
logical :: ok
|
||||
double precision :: phase_Ji, phase_Ik, phase_Ii
|
||||
double precision :: contrib, delta_IJk, HJk, HIk, HIl
|
||||
integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii
|
||||
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
|
||||
integer, allocatable :: idx_sorted_bit(:)
|
||||
integer, external :: get_index_in_psi_det_sorted_bit
|
||||
|
||||
integer :: II, blok
|
||||
|
||||
provide det_cepa0_active delta_cas lambda_mrcc
|
||||
|
||||
if(N_int /= 1) then
|
||||
print *, "mrsc2 experimental N_int==1"
|
||||
stop
|
||||
end if
|
||||
|
||||
i_state = 1
|
||||
delta_sub_ij(:,:,:) = 0d0
|
||||
delta_sub_ii(:,:) = 0d0
|
||||
|
||||
provide mo_bielec_integrals_in_map
|
||||
allocate(idx_sorted_bit(N_det))
|
||||
|
||||
idx_sorted_bit(:) = -1
|
||||
do i=1,N_det_non_ref
|
||||
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
||||
enddo
|
||||
|
||||
!$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_sub_ij, delta_sub_ii)
|
||||
do i=1,N_det_non_ref
|
||||
if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪"
|
||||
do J=1,N_det_ref
|
||||
call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int)
|
||||
if(degree == -1) cycle
|
||||
|
||||
|
||||
do II=1,N_det_ref
|
||||
call apply_excitation(psi_ref(1,1,II),exc_Ji,det_tmp,ok,N_int)
|
||||
!call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,i),exc_Ii,degree,phase_Ii,N_int)
|
||||
|
||||
if(.not. ok) cycle
|
||||
l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
|
||||
if(l == 0) cycle
|
||||
l = idx_sorted_bit(l)
|
||||
|
||||
if(psi_non_ref(1,1,l) /= det_tmp(1,1)) stop "sdf"
|
||||
call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl)
|
||||
|
||||
do k=1,N_det_non_ref
|
||||
call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int)
|
||||
|
||||
det_tmp(:,:) = 0_bit_kind
|
||||
det_tmp2(:,:) = 0_bit_kind
|
||||
|
||||
det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1)))
|
||||
det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1)))
|
||||
ok = (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2))))
|
||||
|
||||
det_tmp(1,1) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,k)), not(active_sorb(2)))
|
||||
det_tmp(1,2) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,i)), not(active_sorb(2)))
|
||||
ok = ok .and. (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2))))
|
||||
|
||||
if(ok) cycle
|
||||
|
||||
|
||||
! call decode_exc(exc_Ii,degree,h1_,p1_,h2_,p2_,s1_,s2_)
|
||||
! call decode_exc(exc_Ik,degree2,h1,p1,h2,p2,s1,s2)
|
||||
!
|
||||
!
|
||||
! det_tmp(:,:) = 0_bit_kind
|
||||
! call set_det_bit(det_tmp, p1, s1)
|
||||
! call set_det_bit(det_tmp, h1, s1)
|
||||
! call set_det_bit(det_tmp, p1_, s1_)
|
||||
! call set_det_bit(det_tmp, h1_, s1_)
|
||||
! if(degree == 2) then
|
||||
! call set_det_bit(det_tmp, p2_, s2_)
|
||||
! call set_det_bit(det_tmp, h2_, s2_)
|
||||
! end if
|
||||
! if(degree2 == 2) then
|
||||
! call set_det_bit(det_tmp, p2, s2)
|
||||
! call set_det_bit(det_tmp, h2, s2)
|
||||
! end if
|
||||
! deg = 0
|
||||
! do ni = 1, N_int
|
||||
! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2))
|
||||
! end do
|
||||
! if(deg == 2*degree2 + 2*degree) cycle
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
! if(degree == -1) cycle
|
||||
call i_h_j(psi_ref(1,1,J), psi_non_ref(1,1,k), N_int, HJk)
|
||||
call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,k), N_int, HIk)
|
||||
if(HJk == 0) cycle
|
||||
!assert HIk == 0
|
||||
delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
|
||||
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
|
||||
if(ok) cycle
|
||||
contrib = delta_IJk * HIl * lambda_mrcc(i_state, l)
|
||||
!$OMP CRITICAL
|
||||
delta_sub_ij(II, i, i_state) += contrib
|
||||
if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
|
||||
delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
|
||||
endif
|
||||
!$OMP END CRITICAL
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!qsd $OMP END PARALLEL DO
|
||||
!$OMP END PARALLEL DO
|
||||
deallocate(idx_sorted_bit)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine set_det_bit(det, p, s)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer(bit_kind),intent(inout) :: det(N_int, 2)
|
||||
integer, intent(in) :: p, s
|
||||
@ -176,6 +358,173 @@ subroutine set_det_bit(det, p, s)
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij_old, (N_det_ref,N_det_non_ref,N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_det_ref,N_states) ]
|
||||
implicit none
|
||||
|
||||
integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni
|
||||
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, x(2), y(2)
|
||||
logical :: ok
|
||||
double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
|
||||
double precision :: contrib
|
||||
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
|
||||
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
|
||||
integer, allocatable :: idx_sorted_bit(:)
|
||||
integer, external :: get_index_in_psi_det_sorted_bit
|
||||
logical, external :: is_in_wavefunction
|
||||
|
||||
delta_ii_old(:,:) = 0
|
||||
delta_ij_old(:,:,:) = 0
|
||||
|
||||
i_state = 1
|
||||
provide mo_bielec_integrals_in_map
|
||||
allocate(idx_sorted_bit(N_det))
|
||||
|
||||
idx_sorted_bit(:) = -1
|
||||
do i=1,N_det_non_ref
|
||||
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
||||
enddo
|
||||
|
||||
!$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij_old, delta_ii_old)
|
||||
do i = 1 , N_det_non_ref
|
||||
if(mod(i,1000) == 0) print *, i, N_det_non_ref
|
||||
do i_I = 1 , N_det_ref
|
||||
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int)
|
||||
if(degree2 == -1) cycle
|
||||
ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
|
||||
|
||||
call decode_exc(exc_iI,degree2,h1,p1,h2,p2,s1,s2)
|
||||
|
||||
call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi)
|
||||
diI = hIi * lambda_mrcc(i_state,i)
|
||||
do J = 1 , N_det_ref !!!
|
||||
call get_excitation(psi_ref(1,1,i_I),psi_ref(1,1,J),exc_IJ,degree,phase_IJ,N_int)
|
||||
call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi)
|
||||
delta_JI = hJi * diI
|
||||
do k = 1 , N_det_non_ref
|
||||
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
|
||||
if(degree == -1) cycle
|
||||
|
||||
call decode_exc(exc_Ik,degree,h1_,p1_,h2_,p2_,s1_,s2_)
|
||||
|
||||
|
||||
det_tmp(:,:) = 0_bit_kind
|
||||
det_tmp2(:,:) = 0_bit_kind
|
||||
|
||||
!!!!!!!!!!!!!!!
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1)))
|
||||
det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1)))
|
||||
ok = (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2))))
|
||||
|
||||
det_tmp(1,1) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,k)), not(active_sorb(2)))
|
||||
det_tmp(1,2) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,i)), not(active_sorb(2)))
|
||||
ok = ok .and. (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2))))
|
||||
if(.not. ok) cycle
|
||||
!if(ok) cycle
|
||||
|
||||
!!!!!!!!!!!!!!
|
||||
|
||||
|
||||
|
||||
! call set_det_bit(det_tmp, p1, s1)
|
||||
!
|
||||
! call set_det_bit(det_tmp, p1_, s1_)
|
||||
!
|
||||
! if(degree == 2) then
|
||||
! call set_det_bit(det_tmp, p2_, s2_)
|
||||
!
|
||||
! end if
|
||||
! if(degree2 == 2) then
|
||||
! call set_det_bit(det_tmp, p2, s2)
|
||||
! end if
|
||||
!
|
||||
! x(:) = 0
|
||||
! do ni=1,N_int
|
||||
! x(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1)))
|
||||
! x(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1)))
|
||||
! end do
|
||||
!
|
||||
!
|
||||
! !det_tmp(:,:) = 0_bit_kind
|
||||
!
|
||||
! call set_det_bit(det_tmp, h1, s1)
|
||||
! call set_det_bit(det_tmp, h1_, s1_)
|
||||
! if(degree == 2) then
|
||||
! call set_det_bit(det_tmp, h2_, s2_)
|
||||
! end if
|
||||
! if(degree2 == 2) then
|
||||
! call set_det_bit(det_tmp, h2, s2)
|
||||
! end if
|
||||
!
|
||||
! y(1) = -x(1)
|
||||
! y(2) = -x(2)
|
||||
! do ni=1,N_int
|
||||
! y(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1)))
|
||||
! y(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1)))
|
||||
! end do
|
||||
!
|
||||
! ! print *, x, y
|
||||
!
|
||||
! if(x(1) * y(1) /= 0) cycle
|
||||
! if(x(2) * y(2) /= 0) cycle
|
||||
!
|
||||
!
|
||||
!
|
||||
! deg = 0
|
||||
! do ni = 1, N_int
|
||||
! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2))
|
||||
! end do
|
||||
! if(deg /= 2*degree2 + 2*degree) cycle
|
||||
|
||||
|
||||
|
||||
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
|
||||
|
||||
call get_excitation(psi_non_ref(1,1,i), det_tmp, exc_Ik, degree, phase_al, N_int)
|
||||
|
||||
|
||||
if(.not. ok) cycle
|
||||
if(is_in_wavefunction(det_tmp, N_int)) cycle
|
||||
|
||||
|
||||
call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp,ok,N_int)
|
||||
if(.not. ok) cycle
|
||||
|
||||
call get_excitation(psi_ref(1,1,J), det_tmp, exc_Ik, degree, phase_Jl, N_int)
|
||||
|
||||
l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
|
||||
if(l == 0) cycle
|
||||
l = idx_sorted_bit(get_index_in_psi_det_sorted_bit(det_tmp, N_int))
|
||||
|
||||
if(l ==-1) cycle
|
||||
|
||||
|
||||
call i_h_j(psi_non_ref(1,1,k), psi_ref(1,1,i_I),N_int,HkI)
|
||||
dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik
|
||||
|
||||
|
||||
!$OMP CRITICAL
|
||||
contrib = dkI(i_state) * delta_JI
|
||||
!erro += abs(dkI(i_state) - psi_non_ref_coef(k,i_state) / psi_ref_coef(1,i_state))
|
||||
delta_ij_old(i_I,l,i_state) += contrib
|
||||
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
|
||||
delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(k,i_state)
|
||||
endif
|
||||
!$OMP END CRITICAL
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
deallocate(idx_sorted_bit)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
subroutine apply_excitation(det, exc, res, ok, Nint)
|
||||
use bitmasks
|
||||
|
||||
@ -55,7 +55,7 @@ subroutine mrcepa0_iterations
|
||||
do i = 1, N_det_ref
|
||||
print*,''
|
||||
print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
|
||||
print*,'delta_ii(i,1) = ',delta_ii(i,1)
|
||||
print*,'delta_ii(i,1) = ',delta_cas(i,i,1)
|
||||
enddo
|
||||
print*,'------------'
|
||||
enddo
|
||||
|
||||
Loading…
Reference in New Issue
Block a user