LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-22 18:57:31 +02:00
Code Issues Releases Wiki Activity

forgotten files

Browse Source
This commit is contained in:
Emmanuel Giner 2016-11-02 16:01:01 +01:00
parent bd91472407
commit 124d918021
33 changed files with 1585 additions and 285 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
plugins/FOBOCI/EZFIO.cfg
View File

@ -19,10 +19,15 @@ default: 0.00001
[do_it_perturbative] [do_it_perturbative]
type: logical type: logical
doc: if true, you do the FOBOCI calculation perturbatively doc: if true, when a given 1h or 1p determinant is not selected because of its perturbation estimate, then if its coefficient is lower than threshold_perturbative, it is acounted in the FOBOCI differential density matrices
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: .False. default: .False.
[threshold_perturbative]
type: double precision
doc: when do_it_perturbative is True, threshold_perturbative select if a given determinant ia selected or not for beign taken into account in the FOBO-SCF treatment. In practive, if the coefficient is larger then threshold_perturbative it means that it not selected as the perturbation should not be too importan. A value of 0.01 is in general OK.
interface: ezfio,provider,ocaml
default: 0.001
[speed_up_convergence_foboscf] [speed_up_convergence_foboscf]
type: logical type: logical
@ -49,3 +54,9 @@ doc: if true, you do all 2p type excitation on the LMCT
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: .True. default: .True.
[selected_fobo_ci]
type: logical
doc: if true, for each CI step you will run a CIPSI calculation that stops at pt2_max
interface: ezfio,provider,ocaml
default: .False.

202
plugins/FOBOCI/SC2_1h1p.irp.f
View File

@ -158,6 +158,7 @@ subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,conve
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma} ! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
! diag_H_elements(index_hf) += total_corr_e_2h2p ! diag_H_elements(index_hf) += total_corr_e_2h2p
return
c_ref = c_ref * c_ref c_ref = c_ref * c_ref
print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf) print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf)
do i = 1, n_singles do i = 1, n_singles
@ -186,6 +187,186 @@ subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,conve
end end
subroutine dressing_1h1p_by_2h2p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: diag_H_elements(dim_in)
double precision, intent(in) :: convergence
integer :: i,j,k,l
integer :: r,s,i0,j0,r0,s0
integer :: n_singles
integer :: index_singles(sze),hole_particles_singles(sze,3)
integer :: n_doubles
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
integer :: index_hf
double precision :: e_corr_singles(mo_tot_num,2)
double precision :: e_corr_doubles(mo_tot_num)
double precision :: e_corr_singles_total(2)
double precision :: e_corr_doubles_1h1p
integer :: exc(0:2,2,2),degree
integer :: h1,h2,p1,p2,s1,s2
integer :: other_spin(2)
double precision :: phase
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i_ok
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz
double precision :: hij,c_ref,contrib
integer :: iorb
other_spin(1) = 2
other_spin(2) = 1
n_singles = 0
n_doubles = 0
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij)
diag_H_elements(i) = hij
if(degree == 0)then
index_hf = i
else if (degree == 1)then
n_singles +=1
index_singles(n_singles) = i
! h1 = inactive orbital of the hole
hole_particles_singles(n_singles,1) = h1
! p1 = virtual orbital of the particle
hole_particles_singles(n_singles,2) = p1
! s1 = spin of the electron excited
hole_particles_singles(n_singles,3) = s1
else if (degree == 2)then
n_doubles +=1
index_doubles(n_doubles) = i
! h1 = inactive orbital of the hole (beta of course)
hole_particles_doubles(n_doubles,1) = h1
! p1 = virtual orbital of the particle (alpha of course)
hole_particles_doubles(n_doubles,2) = p2
else
print*,'PB !! found out other thing than a single or double'
print*,'stopping ..'
stop
endif
enddo
double precision :: delta_e
double precision :: coef_ijrs
diag_H_elements = 0.d0
do i0 = 1, n_core_inact_orb
i= list_core_inact(i0)
do j0 = i0+1, n_core_inact_orb
j = list_core_inact(j0)
print*, i,j
do r0 = 1, n_virt_orb
r = list_virt(r0)
do s0 = r0+1, n_virt_orb
s = list_virt(s0)
!!! alpha (i-->r) / beta (j-->s)
s1 = 1
s2 = 2
key_tmp = ref_bitmask
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
coef_ijrs = hij/delta_e
do k = 1, n_singles
l = index_singles(k)
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
diag_H_elements(l) += coef_ijrs * hij
enddo
!if(i>j.and.r>s)then
!! alpha (i-->r) / alpha (j-->s)
s1 = 1
s2 = 1
key_tmp = ref_bitmask
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
coef_ijrs = hij/delta_e
do k = 1, n_singles
l = index_singles(k)
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
diag_H_elements(l) += coef_ijrs * hij
enddo
!! beta (i-->r) / beta (j-->s)
s1 = 2
s2 = 2
key_tmp = ref_bitmask
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
if(i_ok .ne.1)then
print*, 'pb !!'
stop
endif
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
coef_ijrs = hij/delta_e
do k = 1, n_singles
l = index_singles(k)
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
diag_H_elements(l) += coef_ijrs * hij
enddo
!endif
enddo
enddo
enddo
enddo
c_ref = 1.d0/u_in(index_hf,1)
do k = 1, n_singles
l = index_singles(k)
diag_H_elements(0) -= diag_H_elements(l)
enddo
! do k = 1, n_doubles
! l = index_doubles(k)
! diag_H_elements(0) += diag_H_elements(l)
! enddo
end
subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence) subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence)
use bitmasks use bitmasks
implicit none implicit none
@ -478,11 +659,13 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c
double precision, intent(inout) :: u_in(dim_in,N_st) double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st) double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: diag_H_elements(dim_in) double precision, intent(out) :: diag_H_elements(dim_in)
double precision :: extra_diag_H_elements(dim_in)
double precision, intent(in) :: convergence double precision, intent(in) :: convergence
integer :: i,j,iter integer :: i,j,iter
do iter = 1, 1 do iter = 1, 1
call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) ! call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
if(sze<=N_det_max_jacobi)then call dressing_1h1p_by_2h2p(dets_in,u_in,extra_diag_H_elements,dim_in,sze,N_st,Nint,convergence)
! if(sze<=N_det_max_jacobi)then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
do j=1,sze do j=1,sze
@ -490,9 +673,14 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c
H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
enddo enddo
enddo enddo
do i = 1,sze H_matrix_tmp(1,1) += extra_diag_H_elements(1)
H_matrix_tmp(i,i) = diag_H_elements(i) do i = 2,sze
H_matrix_tmp(1,i) += extra_diag_H_elements(i)
H_matrix_tmp(i,1) += extra_diag_H_elements(i)
enddo enddo
!do i = 1,sze
! H_matrix_tmp(i,i) = diag_H_elements(i)
!enddo
call lapack_diag(eigenvalues,eigenvectors, & call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_tmp,size(H_matrix_all_dets,1),sze) H_matrix_tmp,size(H_matrix_all_dets,1),sze)
do j=1,min(N_states_diag,sze) do j=1,min(N_states_diag,sze)
@ -502,9 +690,9 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c
energies(j) = eigenvalues(j) energies(j) = eigenvalues(j)
enddo enddo
deallocate (H_matrix_tmp, eigenvalues, eigenvectors) deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
else ! else
call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants) ! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants)
endif ! endif
print*,'E = ',energies(1) + nuclear_repulsion print*,'E = ',energies(1) + nuclear_repulsion
enddo enddo

25
plugins/FOBOCI/all_singles.irp.f
View File

@ -1,13 +1,25 @@
subroutine all_single subroutine all_single(e_pt2)
implicit none implicit none
double precision, intent(in) :: e_pt2
integer :: i,k integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree integer :: N_st, degree
double precision,allocatable :: E_before(:) double precision,allocatable :: E_before(:)
N_st = N_states N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
selection_criterion = 0.d0 if(.not.selected_fobo_ci)then
soft_touch selection_criterion selection_criterion = 0.d0
soft_touch selection_criterion
else
selection_criterion = 0.1d0
selection_criterion_factor = 0.01d0
selection_criterion_min = selection_criterion
soft_touch selection_criterion
endif
print*, 'e_pt2 = ',e_pt2
pt2_max = 0.15d0 * e_pt2
soft_touch pt2_max
print*, 'pt2_max = ',pt2_max
threshold_davidson = 1.d-9 threshold_davidson = 1.d-9
soft_touch threshold_davidson davidson_criterion soft_touch threshold_davidson davidson_criterion
i = 0 i = 0
@ -17,6 +29,8 @@ subroutine all_single
print*,'pt2_max = ',pt2_max print*,'pt2_max = ',pt2_max
print*,'N_det_generators = ',N_det_generators print*,'N_det_generators = ',N_det_generators
pt2=-1.d0 pt2=-1.d0
print*, 'ref_bitmask_energy =',ref_bitmask_energy
print*, 'CI_expectation_value =',CI_expectation_value(1)
E_before = ref_bitmask_energy E_before = ref_bitmask_energy
print*,'Initial Step ' print*,'Initial Step '
@ -29,7 +43,7 @@ subroutine all_single
print*,'S^2 = ',CI_eigenvectors_s2(i) print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo enddo
n_det_max = 100000 n_det_max = 100000
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > dabs(pt2_max))
i += 1 i += 1
print*,'-----------------------' print*,'-----------------------'
print*,'i = ',i print*,'i = ',i
@ -39,6 +53,8 @@ subroutine all_single
print*,'E = ',CI_energy(1) print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1) print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1) print*,'E+PT2 = ',E_before + pt2(1)
print*,'pt2_max = ',pt2_max
print*, maxval(abs(pt2(1:N_st))) > dabs(pt2_max)
if(N_states_diag.gt.1)then if(N_states_diag.gt.1)then
print*,'Variational Energy difference' print*,'Variational Energy difference'
do i = 2, N_st do i = 2, N_st
@ -53,7 +69,6 @@ subroutine all_single
endif endif
E_before = CI_energy E_before = CI_energy
!!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO !!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO
exit
enddo enddo
! threshold_davidson = 1.d-8 ! threshold_davidson = 1.d-8
! soft_touch threshold_davidson davidson_criterion ! soft_touch threshold_davidson davidson_criterion

4
plugins/FOBOCI/create_1h_or_1p.irp.f
View File

@ -68,7 +68,9 @@ subroutine create_restart_and_1h(i_hole)
SOFT_TOUCH N_det psi_det psi_coef SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates logical :: found_duplicates
if(n_act_orb.gt.1)then
call remove_duplicates_in_psi_det(found_duplicates) call remove_duplicates_in_psi_det(found_duplicates)
endif
end end
subroutine create_restart_and_1p(i_particle) subroutine create_restart_and_1p(i_particle)
@ -213,6 +215,8 @@ subroutine create_restart_1h_1p(i_hole,i_part)
SOFT_TOUCH N_det psi_det psi_coef SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates logical :: found_duplicates
if(n_act_orb.gt.1)then
call remove_duplicates_in_psi_det(found_duplicates) call remove_duplicates_in_psi_det(found_duplicates)
endif
end end

160
plugins/FOBOCI/dress_simple.irp.f
View File

@ -72,20 +72,21 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen
end end
subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
use bitmasks use bitmasks
implicit none implicit none
double precision, intent(in) :: threshold double precision, intent(in) :: threshold
logical, intent(out) :: is_ok,exit_loop double precision, intent(out):: e_pt2
logical, intent(out) :: is_ok,exit_loop,is_ok_perturbative
logical, intent(in) :: verbose logical, intent(in) :: verbose
integer :: l,k,m integer :: l,k,m
double precision,allocatable :: dressed_H_matrix(:,:) double precision,allocatable :: dressed_H_matrix(:,:)
double precision,allocatable :: psi_coef_diagonalized_tmp(:,:) double precision, allocatable :: psi_coef_diagonalized_tmp(:,:)
integer(bit_kind), allocatable :: psi_det_generators_input(:,:,:) integer(bit_kind), allocatable :: psi_det_generators_input(:,:,:)
double precision :: hij
allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators)) allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators),psi_coef_diagonalized_tmp(N_det_generators,N_states))
allocate(psi_coef_diagonalized_tmp(N_det_generators,N_states))
dressed_H_matrix = 0.d0 dressed_H_matrix = 0.d0
do k = 1, N_det_generators do k = 1, N_det_generators
do l = 1, N_int do l = 1, N_int
@ -94,9 +95,20 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop)
enddo enddo
enddo enddo
!call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input) !call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input)
call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop,is_ok_perturbative)
if(do_it_perturbative)then !do m = 1, N_states
if(is_ok)then ! do k = 1, N_det_generators
! do l = 1, N_int
! psi_selectors(l,1,k) = psi_det_generators_input(l,1,k)
! psi_selectors(l,2,k) = psi_det_generators_input(l,2,k)
! enddo
! psi_selectors_coef(k,m) = psi_coef_diagonalized_tmp(k,m)
! enddo
!enddo
!soft_touch psi_selectors psi_selectors_coef
!if(do_it_perturbative)then
print*, 'is_ok_perturbative',is_ok_perturbative
if(is_ok.or.is_ok_perturbative)then
N_det = N_det_generators N_det = N_det_generators
do m = 1, N_states do m = 1, N_states
do k = 1, N_det_generators do k = 1, N_det_generators
@ -105,11 +117,19 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop)
psi_det(l,2,k) = psi_det_generators_input(l,2,k) psi_det(l,2,k) = psi_det_generators_input(l,2,k)
enddo enddo
psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m) psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m)
print*, 'psi_coef(k,m)',psi_coef(k,m)
enddo
enddo
soft_touch psi_det psi_coef N_det
e_pt2 = 0.d0
do m =1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators_input(1,1,m),psi_det_generators_input(1,1,l),N_int,hij) ! Fill the zeroth order H matrix
e_pt2 += (dressed_H_matrix(m,l) - hij)* psi_coef_diagonalized_tmp(m,1)* psi_coef_diagonalized_tmp(l,1)
enddo enddo
enddo enddo
touch psi_coef psi_det N_det
endif endif
endif !endif
deallocate(psi_det_generators_input,dressed_H_matrix,psi_coef_diagonalized_tmp) deallocate(psi_det_generators_input,dressed_H_matrix,psi_coef_diagonalized_tmp)
@ -118,14 +138,14 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop)
end end
subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop,is_ok_perturbative)
use bitmasks use bitmasks
implicit none implicit none
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
integer, intent(in) :: Ndet_generators integer, intent(in) :: Ndet_generators
double precision, intent(in) :: threshold double precision, intent(in) :: threshold
logical, intent(in) :: verbose logical, intent(in) :: verbose
logical, intent(out) :: is_ok,exit_loop logical, intent(out) :: is_ok,exit_loop,is_ok_perturbative
double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states) double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states)
double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators) double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators)
@ -309,10 +329,124 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
exit exit
endif endif
enddo enddo
if(.not.is_ok)then
is_ok_perturbative = .True.
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
do k = 1, N_states
print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative
if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold_perturbative)then
is_ok_perturbative = .False.
exit
endif
enddo
if(.not.is_ok_perturbative)then
exit
endif
enddo
endif
if(verbose)then if(verbose)then
print*,'is_ok = ',is_ok print*,'is_ok = ',is_ok
print*,'is_ok_perturbative = ',is_ok_perturbative
endif endif
end end
subroutine fill_H_apply_buffer_no_selection_first_order_coef(n_selected,det_buffer,Nint,iproc)
use bitmasks
implicit none
BEGIN_DOC
! Fill the H_apply buffer with determiants for CISD
END_DOC
integer, intent(in) :: n_selected, Nint, iproc
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k
integer :: new_size
PROVIDE H_apply_buffer_allocated
call omp_set_lock(H_apply_buffer_lock(1,iproc))
new_size = H_apply_buffer(iproc)%N_det + n_selected
if (new_size > H_apply_buffer(iproc)%sze) then
call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc)
endif
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
do i=1,n_selected
do j=1,N_int
H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i)
H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i)
enddo
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num)
enddo
double precision :: i_H_psi_array(N_states),h,diag_H_mat_elem_fock,delta_e
do i=1,N_selected
call i_H_psi(det_buffer(1,1,i),psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_psi_array)
call i_H_j(det_buffer(1,1,i),det_buffer(1,1,i),N_int,h)
do j=1,N_states
delta_e = -1.d0 /(h - CI_expectation_value(j))
H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = i_H_psi_array(j) * delta_e
enddo
enddo
H_apply_buffer(iproc)%N_det = new_size
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
call omp_unset_lock(H_apply_buffer_lock(1,iproc))
end
subroutine make_s2_eigenfunction_first_order
implicit none
integer :: i,j,k
integer :: smax, s
integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:)
integer :: N_det_new
integer, parameter :: bufsze = 1000
logical, external :: is_in_wavefunction
allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) )
smax = 1
N_det_new = 0
do i=1,N_occ_pattern
call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int)
s += 1
if (s > smax) then
deallocate(d)
allocate ( d(N_int,2,s) )
smax = s
endif
call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int)
do j=1,s
if (.not. is_in_wavefunction(d(1,1,j), N_int) ) then
N_det_new += 1
do k=1,N_int
det_buffer(k,1,N_det_new) = d(k,1,j)
det_buffer(k,2,N_det_new) = d(k,2,j)
enddo
if (N_det_new == bufsze) then
call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0)
N_det_new = 0
endif
endif
enddo
enddo
if (N_det_new > 0) then
call fill_H_apply_buffer_no_selection_first_order_coef(N_det_new,det_buffer,N_int,0)
call copy_H_apply_buffer_to_wf
SOFT_TOUCH N_det psi_coef psi_det
endif
deallocate(d,det_buffer)
call write_int(output_determinants,N_det_new, 'Added deteminants for S^2')
end

14
plugins/FOBOCI/fobo_scf.irp.f
View File

@ -1,8 +1,13 @@
program foboscf program foboscf
implicit none implicit none
call run_prepare !if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Read" )then
! disk_access_ao_integrals = "Write"
! touch disk_access_ao_integrals
!endif
!print*, 'disk_access_ao_integrals',disk_access_ao_integrals
no_oa_or_av_opt = .True. no_oa_or_av_opt = .True.
touch no_oa_or_av_opt touch no_oa_or_av_opt
call run_prepare
call routine_fobo_scf call routine_fobo_scf
call save_mos call save_mos
@ -10,8 +15,8 @@ end
subroutine run_prepare subroutine run_prepare
implicit none implicit none
no_oa_or_av_opt = .False. ! no_oa_or_av_opt = .False.
touch no_oa_or_av_opt ! touch no_oa_or_av_opt
call damping_SCF call damping_SCF
call diag_inactive_virt_and_update_mos call diag_inactive_virt_and_update_mos
end end
@ -27,6 +32,7 @@ subroutine routine_fobo_scf
print*,'*******************************************************************************' print*,'*******************************************************************************'
print*,'*******************************************************************************' print*,'*******************************************************************************'
print*,'FOBO-SCF Iteration ',i print*,'FOBO-SCF Iteration ',i
print*, 'ao_bielec_integrals_in_map = ',ao_bielec_integrals_in_map
print*,'*******************************************************************************' print*,'*******************************************************************************'
print*,'*******************************************************************************' print*,'*******************************************************************************'
if(speed_up_convergence_foboscf)then if(speed_up_convergence_foboscf)then
@ -46,7 +52,7 @@ subroutine routine_fobo_scf
soft_touch threshold_lmct threshold_mlct soft_touch threshold_lmct threshold_mlct
endif endif
endif endif
call FOBOCI_lmct_mlct_old_thr call FOBOCI_lmct_mlct_old_thr(i)
call save_osoci_natural_mos call save_osoci_natural_mos
call damping_SCF call damping_SCF
call diag_inactive_virt_and_update_mos call diag_inactive_virt_and_update_mos

459
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
View File

@ -1,7 +1,8 @@
subroutine FOBOCI_lmct_mlct_old_thr subroutine FOBOCI_lmct_mlct_old_thr(iter)
use bitmasks use bitmasks
implicit none implicit none
integer, intent(in) :: iter
integer :: i,j,k,l integer :: i,j,k,l
integer(bit_kind),allocatable :: unpaired_bitmask(:,:) integer(bit_kind),allocatable :: unpaired_bitmask(:,:)
integer, allocatable :: occ(:,:) integer, allocatable :: occ(:,:)
@ -10,7 +11,7 @@ subroutine FOBOCI_lmct_mlct_old_thr
logical :: test_sym logical :: test_sym
double precision :: thr,hij double precision :: thr,hij
double precision, allocatable :: dressing_matrix(:,:) double precision, allocatable :: dressing_matrix(:,:)
logical :: verbose,is_ok logical :: verbose,is_ok,is_ok_perturbative
verbose = .True. verbose = .True.
thr = 1.d-12 thr = 1.d-12
allocate(unpaired_bitmask(N_int,2)) allocate(unpaired_bitmask(N_int,2))
@ -46,89 +47,45 @@ subroutine FOBOCI_lmct_mlct_old_thr
i_hole_osoci = list_inact(i) i_hole_osoci = list_inact(i)
print*,'--------------------------' print*,'--------------------------'
! First set the current generators to the one of restart ! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_hole_osoci,thr,test_sym) call check_symetry(i_hole_osoci,thr,test_sym)
if(.not.test_sym)cycle if(.not.test_sym)cycle
call set_generators_to_generators_restart
call set_psi_det_to_generators
print*,'i_hole_osoci = ',i_hole_osoci print*,'i_hole_osoci = ',i_hole_osoci
call create_restart_and_1h(i_hole_osoci) call create_restart_and_1h(i_hole_osoci)
call set_generators_to_psi_det call set_generators_to_psi_det
print*,'Passed set generators' print*,'Passed set generators'
call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask)
call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop) double precision :: e_pt2
call is_a_good_candidate(threshold_lmct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
print*,'is_ok = ',is_ok print*,'is_ok = ',is_ok
if(.not.is_ok)cycle if(is_ok)then
allocate(dressing_matrix(N_det_generators,N_det_generators)) allocate(dressing_matrix(N_det_generators,N_det_generators))
dressing_matrix = 0.d0 dressing_matrix = 0.d0
if(.not.do_it_perturbative)then do k = 1, N_det_generators
do l = 1, N_det_generators
do k = 1, N_det_generators call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
do l = 1, N_det_generators dressing_matrix(k,l) = hkl
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) enddo
dressing_matrix(k,l) = hkl
enddo enddo
enddo hkl = dressing_matrix(1,1)
hkl = dressing_matrix(1,1) do k = 1, N_det_generators
do k = 1, N_det_generators dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
dressing_matrix(k,k) = dressing_matrix(k,k) - hkl enddo
enddo print*,'Naked matrix'
print*,'Naked matrix' do k = 1, N_det_generators
do k = 1, N_det_generators write(*,'(100(F12.5,X))')dressing_matrix(k,:)
write(*,'(100(F12.5,X))')dressing_matrix(k,:) enddo
enddo
! Do all the single excitations on top of the CAS and 1h determinants
! Do all the single excitations on top of the CAS and 1h determinants call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single(e_pt2)
call all_single call make_s2_eigenfunction_first_order
call make_s2_eigenfunction threshold_davidson = 1.d-6
call diagonalize_ci soft_touch threshold_davidson davidson_criterion
! if(dressing_2h2p)then call diagonalize_ci
! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct)
! endif
! ! Change the mask of the holes and particles to perform all the
! ! double excitations that starts from the active space in order
! ! to introduce the Coulomb hole in the active space
! ! These are the 1h2p excitations that have the i_hole_osoci hole in common
! ! and the 2p if there is more than one electron in the active space
! do k = 1, N_int
! zero_bitmask(k,1) = 0_bit_kind
! zero_bitmask(k,2) = 0_bit_kind
! enddo
! ! hole is possible only in the orbital i_hole_osoci
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
! ! and in the active space
! do k = 1, n_act_orb
! call set_bit_to_integer(list_act(k),zero_bitmask(1,1),N_int)
! call set_bit_to_integer(list_act(k),zero_bitmask(1,2),N_int)
! enddo
! call set_bitmask_hole_as_input(zero_bitmask)
! call set_bitmask_particl_as_input(reunion_of_bitmask)
! call all_1h2p
! call diagonalize_CI_SC2
! call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
! ! Change the mask of the holes and particles to perform all the
! ! double excitations that from the orbital i_hole_osoci
! do k = 1, N_int
! zero_bitmask(k,1) = 0_bit_kind
! zero_bitmask(k,2) = 0_bit_kind
! enddo
! ! hole is possible only in the orbital i_hole_osoci
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
! call set_bitmask_hole_as_input(zero_bitmask)
! call set_bitmask_particl_as_input(reunion_of_bitmask)
! call set_psi_det_to_generators
! call all_2h2p
! call diagonalize_CI_SC2
double precision :: hkl double precision :: hkl
call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
hkl = dressing_matrix(1,1) hkl = dressing_matrix(1,1)
@ -139,7 +96,10 @@ subroutine FOBOCI_lmct_mlct_old_thr
do k = 1, N_det_generators do k = 1, N_det_generators
write(*,'(100(F12.5,X))')dressing_matrix(k,:) write(*,'(100(F12.5,X))')dressing_matrix(k,:)
enddo enddo
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) deallocate(dressing_matrix)
else
if(.not.do_it_perturbative)cycle
if(.not. is_ok_perturbative)cycle
endif endif
call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci)
@ -148,7 +108,6 @@ subroutine FOBOCI_lmct_mlct_old_thr
norm_total(k) += norm_tmp(k) norm_total(k) += norm_tmp(k)
enddo enddo
call update_density_matrix_osoci call update_density_matrix_osoci
deallocate(dressing_matrix)
enddo enddo
if(.True.)then if(.True.)then
@ -163,9 +122,9 @@ subroutine FOBOCI_lmct_mlct_old_thr
print*,'--------------------------' print*,'--------------------------'
! First set the current generators to the one of restart ! First set the current generators to the one of restart
call check_symetry(i_particl_osoci,thr,test_sym) call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
call set_generators_to_generators_restart call set_generators_to_generators_restart
call set_psi_det_to_generators call set_psi_det_to_generators
if(.not.test_sym)cycle
print*,'i_particl_osoci= ',i_particl_osoci print*,'i_particl_osoci= ',i_particl_osoci
! Initialize the bitmask to the restart ones ! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones call initialize_bitmask_to_restart_ones
@ -181,32 +140,33 @@ subroutine FOBOCI_lmct_mlct_old_thr
call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask)
!! ! so all the mono excitation on the new generators !! ! so all the mono excitation on the new generators
call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop) call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
print*,'is_ok = ',is_ok print*,'is_ok = ',is_ok
if(.not. is_ok)then if(is_ok)then
allocate(dressing_matrix(N_det_generators,N_det_generators))
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
deallocate(dressing_matrix)
else
if(exit_loop)then if(exit_loop)then
call set_generators_to_generators_restart
call set_psi_det_to_generators
exit exit
else else
cycle if(.not.do_it_perturbative)cycle
if(.not. is_ok_perturbative)cycle
endif endif
endif
allocate(dressing_matrix(N_det_generators,N_det_generators))
if(.not.do_it_perturbative)then
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
call all_single
call make_s2_eigenfunction
call diagonalize_ci
! if(dressing_2h2p)then
! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct)
! endif
endif endif
call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
do k = 1, N_states do k = 1, N_states
@ -214,7 +174,6 @@ subroutine FOBOCI_lmct_mlct_old_thr
norm_total(k) += norm_tmp(k) norm_total(k) += norm_tmp(k)
enddo enddo
call update_density_matrix_osoci call update_density_matrix_osoci
deallocate(dressing_matrix)
enddo enddo
endif endif
@ -376,3 +335,303 @@ subroutine FOBOCI_lmct_old
enddo enddo
print*,'accu = ',accu print*,'accu = ',accu
end end
subroutine FOBOCI_lmct_mlct_old_thr_restart(iter)
use bitmasks
implicit none
integer, intent(in) :: iter
integer :: i,j,k,l
integer(bit_kind),allocatable :: unpaired_bitmask(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha, n_occ_beta
double precision :: norm_tmp(N_states),norm_total(N_states)
logical :: test_sym
double precision :: thr,hij
double precision, allocatable :: dressing_matrix(:,:)
logical :: verbose,is_ok,is_ok_perturbative
verbose = .True.
thr = 1.d-12
allocate(unpaired_bitmask(N_int,2))
allocate (occ(N_int*bit_kind_size,2))
do i = 1, N_int
unpaired_bitmask(i,1) = unpaired_alpha_electrons(i)
unpaired_bitmask(i,2) = unpaired_alpha_electrons(i)
enddo
norm_total = 0.d0
call initialize_density_matrix_osoci
call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int)
print*,''
print*,''
print*,'mulliken spin population analysis'
accu =0.d0
do i = 1, nucl_num
accu += mulliken_spin_densities(i)
print*,i,nucl_charge(i),mulliken_spin_densities(i)
enddo
print*,''
print*,''
print*,'DOING FIRST LMCT !!'
print*,'Threshold_lmct = ',threshold_lmct
integer(bit_kind) , allocatable :: zero_bitmask(:,:)
integer(bit_kind) , allocatable :: psi_singles(:,:,:)
logical :: lmct
double precision, allocatable :: psi_singles_coef(:,:)
logical :: exit_loop
allocate( zero_bitmask(N_int,2) )
if(iter.ne.1)then
do i = 1, n_inact_orb
lmct = .True.
integer :: i_hole_osoci
i_hole_osoci = list_inact(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call check_symetry(i_hole_osoci,thr,test_sym)
if(.not.test_sym)cycle
call set_generators_to_generators_restart
call set_psi_det_to_generators
print*,'i_hole_osoci = ',i_hole_osoci
call create_restart_and_1h(i_hole_osoci)
call set_generators_to_psi_det
print*,'Passed set generators'
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
double precision :: e_pt2
call is_a_good_candidate(threshold_lmct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
print*,'is_ok = ',is_ok
if(is_ok)then
allocate(dressing_matrix(N_det_generators,N_det_generators))
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
hkl = dressing_matrix(1,1)
do k = 1, N_det_generators
dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
enddo
print*,'Naked matrix'
do k = 1, N_det_generators
write(*,'(100(F12.5,X))')dressing_matrix(k,:)
enddo
! Do all the single excitations on top of the CAS and 1h determinants
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
double precision :: hkl
call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
hkl = dressing_matrix(1,1)
do k = 1, N_det_generators
dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
enddo
print*,'Dressed matrix'
do k = 1, N_det_generators
write(*,'(100(F12.5,X))')dressing_matrix(k,:)
enddo
deallocate(dressing_matrix)
else
if(.not.do_it_perturbative)cycle
if(.not. is_ok_perturbative)cycle
endif
call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci)
do k = 1, N_states
print*,'norm_tmp = ',norm_tmp(k)
norm_total(k) += norm_tmp(k)
enddo
call update_density_matrix_osoci
enddo
else
double precision :: array_dm(mo_tot_num)
call read_dm_from_lmct(array_dm)
call update_density_matrix_beta_osoci_read(array_dm)
endif
if(iter.ne.1)then
if(.True.)then
print*,''
print*,'DOING THEN THE MLCT !!'
print*,'Threshold_mlct = ',threshold_mlct
lmct = .False.
do i = 1, n_virt_orb
integer :: i_particl_osoci
i_particl_osoci = list_virt(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
call set_generators_to_generators_restart
call set_psi_det_to_generators
print*,'i_particl_osoci= ',i_particl_osoci
! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call modify_bitmasks_for_particl(i_particl_osoci)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_hole_as_input(unpaired_bitmask)
!!! call all_single_h_core
call create_restart_and_1p(i_particl_osoci)
!!! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
!!! ! so all the mono excitation on the new generators
call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
print*,'is_ok = ',is_ok
if(is_ok)then
allocate(dressing_matrix(N_det_generators,N_det_generators))
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
deallocate(dressing_matrix)
else
if(exit_loop)then
call set_generators_to_generators_restart
call set_psi_det_to_generators
exit
else
if(.not.do_it_perturbative)cycle
if(.not. is_ok_perturbative)cycle
endif
endif
call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
do k = 1, N_states
print*,'norm_tmp = ',norm_tmp(k)
norm_total(k) += norm_tmp(k)
enddo
call update_density_matrix_osoci
enddo
endif
else
integer :: norb
call read_dm_from_mlct(array_dm,norb)
call update_density_matrix_alpha_osoci_read(array_dm)
do i = norb+1, n_virt_orb
i_particl_osoci = list_virt(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
call set_generators_to_generators_restart
call set_psi_det_to_generators
print*,'i_particl_osoci= ',i_particl_osoci
! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call modify_bitmasks_for_particl(i_particl_osoci)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_hole_as_input(unpaired_bitmask)
!!! call all_single_h_core
call create_restart_and_1p(i_particl_osoci)
!!! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
!!! ! so all the mono excitation on the new generators
call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative)
print*,'is_ok = ',is_ok
if(is_ok)then
allocate(dressing_matrix(N_det_generators,N_det_generators))
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
deallocate(dressing_matrix)
else
if(exit_loop)then
call set_generators_to_generators_restart
call set_psi_det_to_generators
exit
else
if(.not.do_it_perturbative)cycle
if(.not. is_ok_perturbative)cycle
endif
endif
call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
do k = 1, N_states
print*,'norm_tmp = ',norm_tmp(k)
norm_total(k) += norm_tmp(k)
enddo
call update_density_matrix_osoci
enddo
endif
print*,'norm_total = ',norm_total
norm_total = norm_generators_restart
norm_total = 1.d0/norm_total
! call rescale_density_matrix_osoci(norm_total)
double precision :: accu
accu = 0.d0
do i = 1, mo_tot_num
accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i)
enddo
print*,'accu = ',accu
end
subroutine read_dm_from_lmct(array)
implicit none
integer :: i,iunit ,getUnitAndOpen
double precision :: stuff
double precision, intent(out) :: array(mo_tot_num)
character*(128) :: input
input=trim("fort.33")
iunit= getUnitAndOpen(input,'r')
print*, iunit
array = 0.d0
do i = 1, n_inact_orb
read(iunit,*) stuff
print*, list_inact(i),stuff
array(list_inact(i)) = stuff
enddo
end
subroutine read_dm_from_mlct(array,norb)
implicit none
integer :: i,iunit ,getUnitAndOpen
double precision :: stuff
double precision, intent(out) :: array(mo_tot_num)
character*(128) :: input
input=trim("fort.35")
iunit= getUnitAndOpen(input,'r')
integer,intent(out) :: norb
read(iunit,*)norb
print*, iunit
input=trim("fort.34")
iunit= getUnitAndOpen(input,'r')
array = 0.d0
print*, 'norb = ',norb
do i = 1, norb
read(iunit,*) stuff
print*, list_virt(i),stuff
array(list_virt(i)) = stuff
enddo
end

2
plugins/FOBOCI/generators_restart_save.irp.f
View File

@ -9,6 +9,7 @@ BEGIN_PROVIDER [ integer, N_det_generators_restart ]
integer :: i integer :: i
integer, save :: ifirst = 0 integer, save :: ifirst = 0
double precision :: norm double precision :: norm
print*, ' Providing N_det_generators_restart'
if(ifirst == 0)then if(ifirst == 0)then
call ezfio_get_determinants_n_det(N_det_generators_restart) call ezfio_get_determinants_n_det(N_det_generators_restart)
ifirst = 1 ifirst = 1
@ -30,6 +31,7 @@ END_PROVIDER
integer :: i, k integer :: i, k
integer, save :: ifirst = 0 integer, save :: ifirst = 0
double precision, allocatable :: psi_coef_read(:,:) double precision, allocatable :: psi_coef_read(:,:)
print*, ' Providing psi_det_generators_restart'
if(ifirst == 0)then if(ifirst == 0)then
call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
do k = 1, N_int do k = 1, N_int

16
plugins/FOBOCI/hcc_1h1p.irp.f
View File

@ -15,12 +15,12 @@ subroutine routine
call diagonalize_CI call diagonalize_CI
call test_hcc call test_hcc
call test_mulliken call test_mulliken
! call SC2_1h1p(psi_det,psi_coef,energies, & call SC2_1h1p(psi_det,psi_coef,energies, &
! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
allocate(H_matrix(N_det,N_det)) ! allocate(H_matrix(N_det,N_det))
call SC2_1h1p_full(psi_det,psi_coef,energies, & ! call SC2_1h1p_full(psi_det,psi_coef,energies, &
H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) ! H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
deallocate(H_matrix) ! deallocate(H_matrix)
integer :: i,j integer :: i,j
double precision :: accu,coef_hf double precision :: accu,coef_hf
! coef_hf = 1.d0/psi_coef(1,1) ! coef_hf = 1.d0/psi_coef(1,1)
@ -34,12 +34,12 @@ end
subroutine pouet subroutine pouet
implicit none implicit none
double precision :: accu,coef_hf double precision :: accu,coef_hf
! provide one_body_dm_mo_alpha one_body_dm_mo_beta provide one_body_dm_mo_alpha one_body_dm_mo_beta
! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int) ! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int)
! touch one_body_dm_mo_alpha one_body_dm_mo_beta ! touch one_body_dm_mo_alpha one_body_dm_mo_beta
call test_hcc call test_hcc
call test_mulliken call test_mulliken
! call save_wavefunction call save_wavefunction
end end

46
plugins/FOBOCI/routines_foboci.irp.f
View File

@ -218,6 +218,44 @@ subroutine update_density_matrix_osoci
enddo enddo
end
subroutine update_density_matrix_beta_osoci_read(array)
implicit none
BEGIN_DOC
! one_body_dm_mo_alpha_osoci += Delta rho alpha
! one_body_dm_mo_beta_osoci += Delta rho beta
END_DOC
integer :: i,j
integer :: iorb,jorb
double precision :: array(mo_tot_num)
do i = 1, mo_tot_num
j = list_act(1)
one_body_dm_mo_beta_osoci(i,j) += array(i)
one_body_dm_mo_beta_osoci(j,i) += array(i)
one_body_dm_mo_beta_osoci(i,i) += array(i) * array(i)
enddo
end
subroutine update_density_matrix_alpha_osoci_read(array)
implicit none
BEGIN_DOC
! one_body_dm_mo_alpha_osoci += Delta rho alpha
! one_body_dm_mo_beta_osoci += Delta rho beta
END_DOC
integer :: i,j
integer :: iorb,jorb
double precision :: array(mo_tot_num)
do i = 1, mo_tot_num
j = list_act(1)
one_body_dm_mo_alpha_osoci(i,j) += array(i)
one_body_dm_mo_alpha_osoci(j,i) += array(i)
one_body_dm_mo_alpha_osoci(i,i) += array(i) * array(i)
enddo
end end
@ -387,14 +425,14 @@ subroutine save_osoci_natural_mos
print*,'ACTIVE ORBITAL ',iorb print*,'ACTIVE ORBITAL ',iorb
do j = 1, n_inact_orb do j = 1, n_inact_orb
jorb = list_inact(j) jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then if(dabs(tmp(iorb,jorb)).gt.0.0001d0)then
print*,'INACTIVE ' print*,'INACTIVE '
print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
endif endif
enddo enddo
do j = 1, n_virt_orb do j = 1, n_virt_orb
jorb = list_virt(j) jorb = list_virt(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then if(dabs(tmp(iorb,jorb)).gt.0.0001d0)then
print*,'VIRT ' print*,'VIRT '
print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
endif endif
@ -412,6 +450,10 @@ subroutine save_osoci_natural_mos
label = "Natural" label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
!if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then
! disk_access_ao_integrals = "Read"
! touch disk_access_ao_integrals
!endif
!soft_touch mo_coef !soft_touch mo_coef
deallocate(tmp,occ) deallocate(tmp,occ)

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

@ -73,21 +73,21 @@
print*, '1h1p = ',accu print*, '1h1p = ',accu
! 1h1p third order ! 1h1p third order
delta_ij_tmp = 0.d0 if(do_third_order_1h1p)then
call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) delta_ij_tmp = 0.d0
!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) 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 do j = 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)
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)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo
enddo enddo
enddo second_order_pt_new_1h1p(i_state) = accu(i_state)
second_order_pt_new_1h1p(i_state) = accu(i_state) enddo
enddo print*, '1h1p(3)',accu
print*, '1h1p(3)',accu endif
! 2h ! 2h
delta_ij_tmp = 0.d0 delta_ij_tmp = 0.d0

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

@ -287,50 +287,46 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
if (n_holes_act == 0 .and. n_particles_act == 1) then if (n_holes_act == 0 .and. n_particles_act == 1) then
! i_particle_act = particles_active_list_spin_traced(1)
! delta_e_act += one_creat_spin_trace(i_particle_act )
ispin = particle_list_practical(1,1) ispin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1) i_particle_act = particle_list_practical(2,1)
call get_excitation_degree(det_1,det_2,degree,N_int) ! call get_excitation_degree(det_1,det_2,degree,N_int)
if(degree == 1)then ! if(degree == 1)then
call get_excitation(det_1,det_2,exc,degree,phase,N_int) ! call get_excitation(det_1,det_2,exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) ! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
i_hole = list_inact_reverse(h1) ! i_hole = list_inact_reverse(h1)
i_part = list_act_reverse(p1) ! i_part = list_act_reverse(p1)
do i_state = 1, N_states ! do i_state = 1, N_states
delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) ! delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state)
enddo ! enddo
else if (degree == 2)then ! else if (degree == 2)then
do i_state = 1, N_states do i_state = 1, N_states
delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state)
enddo enddo
endif ! endif
else if (n_holes_act == 1 .and. n_particles_act == 0) then else if (n_holes_act == 1 .and. n_particles_act == 0) then
! i_hole_act = holes_active_list_spin_traced(1)
! delta_e_act += one_anhil_spin_trace(i_hole_act )
ispin = hole_list_practical(1,1) ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1) i_hole_act = hole_list_practical(2,1)
call get_excitation_degree(det_1,det_2,degree,N_int) ! call get_excitation_degree(det_1,det_2,degree,N_int)
if(degree == 1)then ! if(degree == 1)then
call get_excitation(det_1,det_2,exc,degree,phase,N_int) ! call get_excitation(det_1,det_2,exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) ! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
i_hole = list_act_reverse(h1) ! i_hole = list_act_reverse(h1)
i_part = list_virt_reverse(p1) ! i_part = list_virt_reverse(p1)
do i_state = 1, N_states ! do i_state = 1, N_states
if(isnan(one_creat_virt(i_hole,i_part,i_state)))then ! if(isnan(one_creat_virt(i_hole,i_part,i_state)))then
print*, i_hole,i_part,i_state ! print*, i_hole,i_part,i_state
call debug_det(det_1,N_int) ! call debug_det(det_1,N_int)
call debug_det(det_2,N_int) ! call debug_det(det_2,N_int)
stop ! stop
endif ! endif
delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state) ! delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state)
enddo ! enddo
else if (degree == 2)then ! else if (degree == 2)then
do i_state = 1, N_states do i_state = 1, N_states
delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state)
enddo enddo
endif ! endif
else if (n_holes_act == 1 .and. n_particles_act == 1) then else if (n_holes_act == 1 .and. n_particles_act == 1) then
! i_hole_act = holes_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
@ -460,21 +456,12 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
double precision :: phase double precision :: phase
call get_excitation_degree(det_1,det_2,degree,N_int) call get_excitation_degree(det_1,det_2,degree,N_int)
if(degree == 1)then if(degree == 1)then
! call debug_det(det_1,N_int)
call get_excitation(det_1,det_2,exc,degree,phase,N_int) call get_excitation(det_1,det_2,exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
i_hole = list_inact_reverse(h1) i_hole = list_inact_reverse(h1)
i_part = list_virt_reverse(p1) i_part = list_virt_reverse(p1)
do i_state = 1, N_states do i_state = 1, N_states
! if(one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1).gt.1.d-10)then ! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state)
! print*, hij, one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) &
! * coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
! print*, coef_array(i_state)* hij*coef_array(i_state)* hij,one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1), &
! coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
! else
delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state)
! endif
enddo enddo
endif endif

2
plugins/Perturbation/EZFIO.cfg
View File

@ -16,4 +16,4 @@ type: Normalized_float
doc: The selection process stops when the energy ratio variational/(variational+PT2) doc: The selection process stops when the energy ratio variational/(variational+PT2)
is equal to var_pt2_ratio is equal to var_pt2_ratio
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 0.75 default: 0.75

54
plugins/Perturbation/pt2_equations.irp.f
View File

@ -142,6 +142,60 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
end end
subroutine pt2_epstein_nesbet_2x2_no_ci_diag($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various N_st states.
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem_fock,delta_e, h
double precision :: i_H_psi_array(N_st)
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
PROVIDE CI_electronic_energy
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
!call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if (i_H_psi_array(i) /= 0.d0) then
delta_e = h - CI_expectation_value(i)
if (delta_e > 0.d0) then
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
else
e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
endif
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
else
e_2_pert(i) = 0.d0
c_pert(i) = 0.d0
H_pert_diag(i) = 0.d0
endif
enddo
end
subroutine pt2_moller_plesset ($arguments) subroutine pt2_moller_plesset ($arguments)
use bitmasks use bitmasks
implicit none implicit none

5
plugins/Properties/hyperfine_constants.irp.f
View File

@ -121,6 +121,11 @@ END_PROVIDER
conversion_factor_mhz_hcc(8) = -606.1958551736545d0 conversion_factor_mhz_hcc(8) = -606.1958551736545d0
conversion_factor_gauss_hcc(8) = -216.30574771560407d0 conversion_factor_gauss_hcc(8) = -216.30574771560407d0
conversion_factor_cm_1_hcc(8) = -202.20517197179822d0 conversion_factor_cm_1_hcc(8) = -202.20517197179822d0
! Phosphore
conversion_factor_mhz_hcc(15) = 1811.0967763744873d0
conversion_factor_gauss_hcc(15) = 646.2445276897648d0
conversion_factor_cm_1_hcc(15) = 604.1170297381395d0
END_PROVIDER END_PROVIDER

1
plugins/Selectors_no_sorted/selectors.irp.f
View File

@ -40,6 +40,7 @@ END_PROVIDER
do k=1,N_states do k=1,N_states
do i=1,N_det_selectors do i=1,N_det_selectors
psi_selectors_coef(i,k) = psi_coef(i,k) psi_selectors_coef(i,k) = psi_coef(i,k)
! print*, 'psi_selectors_coef(i,k) == ',psi_selectors_coef(i,k)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER

110
plugins/loc_cele/loc_exchange_int.irp.f Normal file
View File

@ -0,0 +1,110 @@
program loc_int
implicit none
integer :: i,j,k,l,iorb,jorb
double precision :: exchange_int(mo_tot_num)
integer :: iorder(mo_tot_num)
integer :: indices(mo_tot_num,2)
logical :: list_core_inact_check(mo_tot_num)
integer :: n_rot
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_core_inact_orb
iorb = list_core_inact(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_core_inact_orb
jorb = list_core_inact(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_act_orb
iorb = list_act(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_act_orb
jorb = list_act(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_virt_orb
iorb = list_virt(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_virt_orb
jorb = list_virt(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
call save_mos
end

45
plugins/loc_cele/loc_exchange_int_act.irp.f Normal file
View File

@ -0,0 +1,45 @@
program loc_int
implicit none
integer :: i,j,k,l,iorb,jorb
double precision :: exchange_int(mo_tot_num)
integer :: iorder(mo_tot_num)
integer :: indices(mo_tot_num,2)
logical :: list_core_inact_check(mo_tot_num)
integer :: n_rot
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_act_orb
iorb = list_act(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_act_orb
jorb = list_act(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
call save_mos
end

45
plugins/loc_cele/loc_exchange_int_inact.irp.f Normal file
View File

@ -0,0 +1,45 @@
program loc_int
implicit none
integer :: i,j,k,l,iorb,jorb
double precision :: exchange_int(mo_tot_num)
integer :: iorder(mo_tot_num)
integer :: indices(mo_tot_num,2)
logical :: list_core_inact_check(mo_tot_num)
integer :: n_rot
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_core_inact_orb
iorb = list_core_inact(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_core_inact_orb
jorb = list_core_inact(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
call save_mos
end

47
plugins/loc_cele/loc_exchange_int_virt.irp.f Normal file
View File

@ -0,0 +1,47 @@
program loc_int
implicit none
integer :: i,j,k,l,iorb,jorb
double precision :: exchange_int(mo_tot_num)
integer :: iorder(mo_tot_num)
integer :: indices(mo_tot_num,2)
logical :: list_core_inact_check(mo_tot_num)
integer :: n_rot
indices = 0
list_core_inact_check = .True.
n_rot = 0
do i = 1, n_virt_orb
iorb = list_virt(i)
exchange_int = 0.d0
iorder = 0
print*,''
if(list_core_inact_check(iorb) == .False.)cycle
do j = i+1, n_virt_orb
jorb = list_virt(j)
iorder(jorb) = jorb
exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
enddo
n_rot += 1
call dsort(exchange_int,iorder,mo_tot_num)
indices(n_rot,1) = iorb
indices(n_rot,2) = iorder(1)
list_core_inact_check(iorder(1)) = .False.
print*,indices(n_rot,1),indices(n_rot,2)
print*,''
print*,''
enddo
print*,'****************************'
print*,'-+++++++++++++++++++++++++'
do i = 1, n_rot
iorb = indices(i,1)
jorb = indices(i,2)
print*,iorb,jorb
call mix_mo_jk(iorb,jorb)
enddo
call save_mos
end

13
src/Determinants/EZFIO.cfg
View File

@ -124,3 +124,16 @@ interface: ezfio,provider,ocaml
doc: Energy that should be obtained when truncating the wave function (optional) doc: Energy that should be obtained when truncating the wave function (optional)
type: Energy type: Energy
default: 0. default: 0.
[store_full_H_mat]
type: logical
doc: If True, the Davidson diagonalization is performed by storring the full H matrix up to n_det_max_stored. Be carefull, it can cost a lot of memory but can also save a lot of CPU time
interface: ezfio,provider,ocaml
default: False
[n_det_max_stored]
type: Det_number_max
doc: Maximum number of determinants for which the full H matrix is stored. Be carefull, the memory requested scales as 10*n_det_max_stored**2. For instance, 90000 determinants represent a matrix of size 60 Gb.
interface: ezfio,provider,ocaml
default: 90000

1
src/Determinants/H_apply.irp.f
View File

@ -306,7 +306,6 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
call omp_unset_lock(H_apply_buffer_lock(1,iproc)) call omp_unset_lock(H_apply_buffer_lock(1,iproc))
end end
subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id)
use f77_zmq use f77_zmq
implicit none implicit none

9
src/Determinants/davidson.irp.f
View File

@ -334,6 +334,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
double precision :: to_print(2,N_st) double precision :: to_print(2,N_st)
double precision :: cpu, wall double precision :: cpu, wall
if(store_full_H_mat.and.sze.le.n_det_max_stored)then
provide H_matrix_all_dets
endif
call write_time(iunit) call write_time(iunit)
@ -439,7 +442,11 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
! ---------------------- ! ----------------------
do k=1,N_st do k=1,N_st
call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) if(store_full_H_mat.and.sze.le.n_det_max_stored)then
call H_u_0_stored(W(1,k,iter),U(1,k,iter),H_matrix_all_dets,sze)
else
call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
endif
enddo enddo

25
src/Determinants/diagonalize_CI.irp.f
View File

@ -33,6 +33,14 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_expectation_value, (N_states_diag) ]
implicit none
integer :: i
do i = 1, N_states
call u0_H_u_0(CI_expectation_value(i),psi_coef(1,i),n_det,psi_det,N_int)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ] BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ]
@ -69,10 +77,14 @@ END_PROVIDER
if (diag_algorithm == "Davidson") then if (diag_algorithm == "Davidson") then
print*, '------------- In Davidson '
call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, & call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, &
size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_determinants) size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_determinants)
print*, '------------- Out Davidson '
do j=1,N_states_diag do j=1,N_states_diag
print*, '------------- In S^2'
call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),CI_eigenvectors_s2(j)) call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),CI_eigenvectors_s2(j))
print*, '------------- Out S^2'
enddo enddo
@ -233,16 +245,20 @@ END_PROVIDER
else else
! Sorting the N_states_diag by energy, whatever the S^2 value is !! Sorting the N_states_diag by energy, whatever the S^2 value is
allocate(e_array(n_states_diag),iorder(n_states_diag)) allocate(e_array(n_states_diag),iorder(n_states_diag))
do j = 1, N_states_diag do j = 2, N_states_diag
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) if(store_full_H_mat.and.n_det.le.n_det_max_stored)then
call u_0_H_u_0_stored(e_0,CI_eigenvectors(1,j),H_matrix_all_dets,n_det)
else
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
endif
e_array(j) = e_0 e_array(j) = e_0
iorder(j) = j iorder(j) = j
enddo enddo
call dsort(e_array,iorder,n_states_diag) call dsort(e_array,iorder,n_states_diag)
do j = 1, N_states_diag do j = 2, N_states_diag
CI_electronic_energy(j) = e_array(j) CI_electronic_energy(j) = e_array(j)
do i = 1, N_det do i = 1, N_det
CI_eigenvectors(i,j) = psi_coef(i,iorder(j)) CI_eigenvectors(i,j) = psi_coef(i,iorder(j))
@ -253,6 +269,7 @@ END_PROVIDER
endif endif
deallocate(s2_eigvalues) deallocate(s2_eigvalues)
endif endif
print*, 'out provider'
END_PROVIDER END_PROVIDER

43
src/Determinants/occ_pattern.irp.f
View File

@ -256,27 +256,6 @@ subroutine make_s2_eigenfunction
integer :: N_det_new integer :: N_det_new
integer, parameter :: bufsze = 1000 integer, parameter :: bufsze = 1000
logical, external :: is_in_wavefunction logical, external :: is_in_wavefunction
! return
! !TODO DEBUG
! do i=1,N_det
! do j=i+1,N_det
! s = 0
! do k=1,N_int
! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. &
! (psi_det(k,2,j) /= psi_det(k,2,i))) then
! s=1
! exit
! endif
! enddo
! if ( s == 0 ) then
! print *, 'Error0: det ', j, 'already in wf'
! call debug_det(psi_det(1,1,j),N_int)
! stop
! endif
! enddo
! enddo
! !TODO DEBUG
allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) ) allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) )
smax = 1 smax = 1
@ -308,33 +287,15 @@ subroutine make_s2_eigenfunction
if (N_det_new > 0) then if (N_det_new > 0) then
call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0) call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0)
! call fill_H_apply_buffer_no_selection_first_order_coef(N_det_new,det_buffer,N_int,0)
call copy_H_apply_buffer_to_wf call copy_H_apply_buffer_to_wf
SOFT_TOUCH N_det psi_coef psi_det SOFT_TOUCH N_det psi_coef psi_det
endif endif
deallocate(d,det_buffer) deallocate(d,det_buffer)
! !TODO DEBUG
! do i=1,N_det
! do j=i+1,N_det
! s = 0
! do k=1,N_int
! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. &
! (psi_det(k,2,j) /= psi_det(k,2,i))) then
! s=1
! exit
! endif
! enddo
! if ( s == 0 ) then
! print *, 'Error : det ', j, 'already in wf at ', i
! call debug_det(psi_det(1,1,j),N_int)
! stop
! endif
! enddo
! enddo
! !TODO DEBUG
call write_int(output_determinants,N_det_new, 'Added deteminants for S^2') call write_int(output_determinants,N_det_new, 'Added deteminants for S^2')
end end

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

@ -431,7 +431,7 @@ end
subroutine i_H_j(key_i,key_j,Nint,hij) subroutine i_H_j_new(key_i,key_j,Nint,hij)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -463,6 +463,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
hij = 0.d0 hij = 0.d0
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint) call get_excitation_degree(key_i,key_j,degree,Nint)
integer :: spin
select case (degree) select case (degree)
case (2) case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint) call get_double_excitation(key_i,key_j,exc,phase,Nint)
@ -507,6 +508,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
! Mono alpha ! Mono alpha
m = exc(1,1,1) m = exc(1,1,1)
p = exc(1,2,1) p = exc(1,2,1)
spin = 1
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
@ -534,6 +536,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
! Mono beta ! Mono beta
m = exc(1,1,2) m = exc(1,1,2)
p = exc(1,2,2) p = exc(1,2,2)
spin = 2
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
@ -559,6 +563,154 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
endif endif
hij = phase*(hij + mo_mono_elec_integral(m,p)) hij = phase*(hij + mo_mono_elec_integral(m,p))
case (0)
hij = diag_H_mat_elem(key_i,Nint)
end select
end
subroutine i_H_j(key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
integer :: exc(0:2,2,2)
integer :: degree
double precision :: get_mo_bielec_integral_schwartz
integer :: m,n,p,q
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
!DIR$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
integer :: spin
select case (degree)
case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
if(exc(1,1,1) == exc(1,2,2) )then
hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) ==exc(1,1,2))then
hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
else
hij = phase*get_mo_bielec_integral_schwartz( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map)
endif
else if (exc(0,1,1) == 2) then
! Double alpha
hij = phase*(get_mo_bielec_integral_schwartz( &
exc(1,1,1), &
exc(2,1,1), &
exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - &
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
! Double beta
hij = phase*(get_mo_bielec_integral_schwartz( &
exc(1,1,2), &
exc(2,1,2), &
exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( &
exc(1,1,2), &
exc(2,1,2), &
exc(2,2,2), &
exc(1,2,2) ,mo_integrals_map) )
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
spin = 1
! do k = 1, elec_alpha_num
! i = occ(k,1)
! if (.not.has_mipi(i)) then
! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
! miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
! has_mipi(i) = .True.
! endif
! enddo
! do k = 1, elec_beta_num
! i = occ(k,2)
! if (.not.has_mipi(i)) then
! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
! has_mipi(i) = .True.
! endif
! enddo
!
! do k = 1, elec_alpha_num
! hij = hij + mipi(occ(k,1)) - miip(occ(k,1))
! enddo
! do k = 1, elec_beta_num
! hij = hij + mipi(occ(k,2))
! enddo
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
spin = 2
! do k = 1, elec_beta_num
! i = occ(k,2)
! if (.not.has_mipi(i)) then
! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
! miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
! has_mipi(i) = .True.
! endif
! enddo
! do k = 1, elec_alpha_num
! i = occ(k,1)
! if (.not.has_mipi(i)) then
! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
! has_mipi(i) = .True.
! endif
! enddo
!
! do k = 1, elec_alpha_num
! hij = hij + mipi(occ(k,1))
! enddo
! do k = 1, elec_beta_num
! hij = hij + mipi(occ(k,2)) - miip(occ(k,2))
! enddo
endif
! hij = phase*(hij + mo_mono_elec_integral(m,p))
call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
case (0) case (0)
hij = diag_H_mat_elem(key_i,Nint) hij = diag_H_mat_elem(key_i,Nint)
@ -2182,3 +2334,43 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
deallocate (shortcut, sort_idx, sorted, version) deallocate (shortcut, sort_idx, sorted, version)
end end
subroutine H_u_0_stored(v_0,u_0,hmatrix,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! uses the big_matrix_stored array
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: hmatrix(sze,sze)
double precision, intent(out) :: v_0(sze)
double precision, intent(in) :: u_0(sze)
v_0 = 0.d0
call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
end
subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>
!
! n : number of determinants
!
! uses the big_matrix_stored array
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: hmatrix(sze,sze)
double precision, intent(out) :: e_0
double precision, intent(in) :: u_0(sze)
double precision :: v_0(sze)
double precision :: u_dot_v
e_0 = 0.d0
v_0 = 0.d0
call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
e_0 = u_dot_v(v_0,u_0,sze)
end

42
src/Determinants/utils.irp.f
View File

@ -1,15 +1,21 @@
BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
use bitmasks
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! H matrix on the basis of the slater determinants defined by psi_det ! H matrix on the basis of the slater determinants defined by psi_det
END_DOC END_DOC
integer :: i,j integer :: i,j,k
double precision :: hij double precision :: hij
integer :: degree(N_det),idx(0:N_det)
call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
!$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij) & !$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij,degree,idx,k) &
!$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets) !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets)
do i =1,N_det do i =1,N_det
do j =i,N_det ! call get_excitation_degree_vector(psi_det,psi_det(1,1,i),degree,N_int,N_det,idx)
! do k =1, idx(0)
! j = idx(k)
! if(j.lt.i)cycle
do j = i, N_det
call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
H_matrix_all_dets(i,j) = hij H_matrix_all_dets(i,j) = hij
H_matrix_all_dets(j,i) = hij H_matrix_all_dets(j,i) = hij
@ -18,3 +24,33 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
END_PROVIDER END_PROVIDER
subroutine provide_big_matrix_stored_with_current_dets(sze,dets_in,big_matrix_stored)
use bitmasks
integer, intent(in) :: sze
integer(bit_kind), intent(in) :: dets_in(N_int,2,sze)
double precision, intent(out) :: big_matrix_stored(sze,sze)
integer :: i,j,k
double precision :: hij
integer :: degree(N_det),idx(0:N_det)
call i_H_j(dets_in(1,1,1),dets_in(1,1,1),N_int,hij)
print*, 'providing big_matrix_stored'
print*, n_det_max_stored
!$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij,degree,idx,k) &
!$OMP SHARED (sze, dets_in, N_int,big_matrix_stored)
do i =1,sze
! call get_excitation_degree_vector(dets_in,dets_in(1,1,i),degree,N_int,sze,idx)
! do k =1, idx(0)
! j = idx(k)
do j = i, sze
if(j.lt.i)cycle
call i_H_j(dets_in(1,1,i),dets_in(1,1,j),N_int,hij)
big_matrix_stored(i,j) = hij
big_matrix_stored(j,i) = hij
enddo
enddo
!$OMP END PARALLEL DO
print*, 'big_matrix_stored provided !!'
end

23
src/Integrals_Bielec/EZFIO.cfg
View File

@ -12,6 +12,22 @@ interface: ezfio,provider,ocaml
default: False default: False
ezfio_name: no_vvvv_integrals ezfio_name: no_vvvv_integrals
[write_ao_map_after_transfo]
type: logical
doc: If True, you dump all the ao integrals after having transformed the mo integrals
interface: ezfio,provider,ocaml
default: False
ezfio_name: write_ao_map_after_transfo
[clear_ao_map_after_mo_transfo]
type: logical
doc: If True, you clear all the ao integrals after having done the transformation
interface: ezfio,provider,ocaml
default: False
ezfio_name: clear_ao_map_after_mo_transfo
[no_ivvv_integrals] [no_ivvv_integrals]
type: logical type: logical
doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual index and 1 belonging to the core inactive active orbitals doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual index and 1 belonging to the core inactive active orbitals
@ -19,6 +35,13 @@ interface: ezfio,provider,ocaml
default: False default: False
ezfio_name: no_ivvv_integrals ezfio_name: no_ivvv_integrals
[no_vvv_integrals]
type: logical
doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual orbitals
interface: ezfio,provider,ocaml
default: False
ezfio_name: no_vvv_integrals
[disk_access_mo_integrals] [disk_access_mo_integrals]
type: Disk_access type: Disk_access
doc: Read/Write MO integrals from/to disk [ Write | Read | None ] doc: Read/Write MO integrals from/to disk [ Write | Read | None ]

2
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -349,6 +349,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
integral = ao_bielec_integral(1,1,1,1) integral = ao_bielec_integral(1,1,1,1)
real :: map_mb real :: map_mb
print*, 'read_ao_integrals',read_ao_integrals
print*, 'disk_access_ao_integrals',disk_access_ao_integrals
if (read_ao_integrals) then if (read_ao_integrals) then
integer :: load_ao_integrals integer :: load_ao_integrals
print*,'Reading the AO integrals' print*,'Reading the AO integrals'

22
src/Integrals_Bielec/integrals_3_index.irp.f Normal file
View File

@ -0,0 +1,22 @@
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num_align,mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num_align,mo_tot_num, mo_tot_num)]
implicit none
integer :: i,j,k,l
double precision :: get_mo_bielec_integral_schwartz
double precision :: integral
do k = 1, mo_tot_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
l = j
integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map)
big_array_coulomb_integrals(j,i,k) = integral
l = j
integral = get_mo_bielec_integral_schwartz(i,j,l,k,mo_integrals_map)
big_array_exchange_integrals(j,i,k) = integral
enddo
enddo
enddo
END_PROVIDER

67
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -414,6 +414,73 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map)
deallocate(pairs,hash,iorder,tmp_val) deallocate(pairs,hash,iorder,tmp_val)
end end
subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ki|li>
! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1)
! for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_in_map
integer :: kk
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(k,i,l,i,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
subroutine get_mo_bielec_integrals_exch_ii(k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ki|il>
! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1)
! for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_in_map
integer :: kk
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(k,i,i,l,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
integer*8 function get_mo_map_size() integer*8 function get_mo_map_size()
implicit none implicit none
BEGIN_DOC BEGIN_DOC

65
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -20,6 +20,7 @@ end
BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
use map_module
implicit none implicit none
integer(bit_kind) :: mask_ijkl(N_int,4) integer(bit_kind) :: mask_ijkl(N_int,4)
integer(bit_kind) :: mask_ijk(N_int,3) integer(bit_kind) :: mask_ijk(N_int,3)
@ -40,7 +41,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
if(no_vvvv_integrals)then if(no_vvvv_integrals)then
integer :: i,j,k,l integer :: i,j,k,l
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!!
! (core+inact+act) ^ 4 ! (core+inact+act) ^ 4
! <ii|ii> ! <ii|ii>
print*, '' print*, ''
@ -52,7 +53,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1)
enddo enddo
call add_integrals_to_map(mask_ijkl) call add_integrals_to_map(mask_ijkl)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!!
! (core+inact+act) ^ 2 (virt) ^2 ! (core+inact+act) ^ 2 (virt) ^2
! <iv|iv> = J_iv ! <iv|iv> = J_iv
print*, '' print*, ''
@ -76,17 +77,19 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
mask_ijkl(i,4) = virt_bitmask(i,1) mask_ijkl(i,4) = virt_bitmask(i,1)
enddo enddo
call add_integrals_to_map(mask_ijkl) call add_integrals_to_map(mask_ijkl)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!!
if(.not.no_vvv_integrals)then
print*, '' print*, ''
print*, '<vr|vs>' print*, '<rv|sv> and <rv|vs>'
do i = 1,N_int do i = 1,N_int
mask_ijk(i,1) = virt_bitmask(i,1) mask_ijk(i,1) = virt_bitmask(i,1)
mask_ijk(i,2) = virt_bitmask(i,1) mask_ijk(i,2) = virt_bitmask(i,1)
mask_ijk(i,3) = virt_bitmask(i,1) mask_ijk(i,3) = virt_bitmask(i,1)
enddo enddo
call add_integrals_to_map_three_indices(mask_ijk) call add_integrals_to_map_three_indices(mask_ijk)
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!!
! (core+inact+act) ^ 3 (virt) ^1 ! (core+inact+act) ^ 3 (virt) ^1
! <iv|ii> ! <iv|ii>
print*, '' print*, ''
@ -101,9 +104,9 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!!
! (core+inact+act) ^ 1 (virt) ^3 ! (core+inact+act) ^ 1 (virt) ^3
! <iv|vv> ! <iv|vv>
print*, ''
print*, '<iv|vv>'
if(.not.no_ivvv_integrals)then if(.not.no_ivvv_integrals)then
print*, ''
print*, '<iv|vv>'
do i = 1,N_int do i = 1,N_int
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
mask_ijkl(i,2) = virt_bitmask(i,1) mask_ijkl(i,2) = virt_bitmask(i,1)
@ -116,6 +119,21 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
else else
call add_integrals_to_map(full_ijkl_bitmask_4) call add_integrals_to_map(full_ijkl_bitmask_4)
endif endif
if(write_ao_map_after_transfo)then
call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin')
disk_access_ao_integrals = "Read"
touch disk_access_ao_integrals
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
endif
if(clear_ao_map_after_mo_transfo)then
call clear_ao_map
integer (map_size_kind) :: get_ao_map_size
print*, '^^^^^^^^^^^^^^^^^^^^^'
print *, 'get_ao_map_size',get_ao_map_size
print*, '^^^^^^^^^^^^^^^^^^^^^'
FREE ao_bielec_integrals_in_map
endif
END_PROVIDER END_PROVIDER
subroutine set_integrals_jj_into_map subroutine set_integrals_jj_into_map
@ -391,64 +409,41 @@ subroutine add_integrals_to_map(mask_ijkl)
endif endif
j1 = ishft((l*l-l),-1) j1 = ishft((l*l-l),-1)
do j0 = 1, n_j do j0 = 1, n_j
! print*, 'l :: j0',l
j = list_ijkl(j0,2) j = list_ijkl(j0,2)
! print*, 'j :: 2',j
if (j > l) then if (j > l) then
! print*, 'j>l'
! print*, j,l
exit exit
endif endif
j1 += 1 j1 += 1
do k0 = 1, n_k do k0 = 1, n_k
k = list_ijkl(k0,3) k = list_ijkl(k0,3)
! print*, 'l :: k0',l
! print*, 'k :: 3',j
i1 = ishft((k*k-k),-1) i1 = ishft((k*k-k),-1)
if (i1<=j1) then if (i1<=j1) then
continue continue
else else
! print*, 'k>l'
! print*, k,l
exit exit
endif endif
bielec_tmp_1 = 0.d0 bielec_tmp_1 = 0.d0
do i0 = 1, n_i do i0 = 1, n_i
i = list_ijkl(i0,1) i = list_ijkl(i0,1)
! print*, 'l :: i0',l
! print*, 'i :: 1',i
if (i>k) then if (i>k) then
! print*, 'i>k'
! print*, i,k
exit exit
endif endif
bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0)
! i1+=1 ! i1+=1
enddo enddo
! do i = 1, min(k,j1-i1+list_ijkl(1,1))
! do i = 1, min(k,j1-i1+list_ijkl(1,1)-1)
do i0 = 1, n_i do i0 = 1, n_i
i = list_ijkl(i0,1) i = list_ijkl(i0,1)
if(i> min(k,j1-i1+list_ijkl(1,1)-1))then if(i> min(k,j1-i1+list_ijkl(1,1)-1))then
! if (i>k) then !min(k,j1-i1)
exit exit
endif endif
! print*, i,j,k,l
! print*, k,j1,i1,j1-i1
if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then
cycle cycle
endif endif
! print*, i,j,k,l
n_integrals += 1 n_integrals += 1
buffer_value(n_integrals) = bielec_tmp_1(i) buffer_value(n_integrals) = bielec_tmp_1(i)
!DEC$ FORCEINLINE !DEC$ FORCEINLINE
call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
! if(i==12.and.k==12 .and.j==12.and.l==12)then
! print*, i,j,k,l,buffer_i(n_integrals)
! accu_bis += buffer_value(n_integrals)
! print*, buffer_value(n_integrals),accu_bis
! endif
if (n_integrals == size_buffer) then if (n_integrals == size_buffer) then
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind)) real(mo_integrals_threshold,integral_kind))
@ -631,7 +626,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk)
bielec_tmp_2 = 0.d0 bielec_tmp_2 = 0.d0
do j1 = 1,ao_num do j1 = 1,ao_num
call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
enddo enddo
do j1 = 1,ao_num do j1 = 1,ao_num
kmax = 0 kmax = 0
@ -732,9 +726,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk)
j = list_ijkl(j0,2) j = list_ijkl(j0,2)
do i0 = 1, n_i do i0 = 1, n_i
i = list_ijkl(i0,1) i = list_ijkl(i0,1)
! if(m==2)then
! if(i==j .and. j == k)cycle
! endif
if (i>k) then if (i>k) then
exit exit
endif endif

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

@ -561,3 +561,18 @@ end
subroutine matrix_vector_product(u0,u1,matrix,sze,lda)
implicit none
BEGIN_DOC
! performs u1 += u0 * matrix
END_DOC
integer, intent(in) :: sze,lda
double precision, intent(in) :: u0(sze)
double precision, intent(inout) :: u1(sze)
double precision, intent(in) :: matrix(lda,sze)
integer :: i,j
integer :: incx,incy
incx = 1
incy = 1
call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy)
end