LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-07 14:03:37 +01:00
Code Releases Activity

Merge pull request #4 from kgasperich/features_green

Browse Source 
fixed transformation
This commit is contained in:
Kevin Gasperich 2020-06-11 13:42:50 -05:00 committed by GitHub
commit 35773160b0
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
18 changed files with 3391 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/bitmask/core_inact_act_virt.irp.f
View File

@ -448,7 +448,7 @@ BEGIN_PROVIDER [ integer, n_core_orb_kpts, (kpt_num)]
do k=1,kpt_num
n_core_orb_kpts(k) = 0
kshift = (1-k)*mo_num_per_kpt
kshift = (k-1)*mo_num_per_kpt
do i = 1, mo_num_per_kpt
if(mo_class(i+kshift) == 'Core')then
n_core_orb_kpts(k) += 1
@ -469,7 +469,7 @@ BEGIN_PROVIDER [ integer, n_inact_orb_kpts, (kpt_num)]
do k=1,kpt_num
n_inact_orb_kpts(k) = 0
kshift = (1-k)*mo_num_per_kpt
kshift = (k-1)*mo_num_per_kpt
do i = 1, mo_num_per_kpt
if(mo_class(i+kshift) == 'Inactive')then
n_inact_orb_kpts(k) += 1
@ -490,7 +490,7 @@ BEGIN_PROVIDER [ integer, n_act_orb_kpts, (kpt_num)]
do k=1,kpt_num
n_act_orb_kpts(k) = 0
kshift = (1-k)*mo_num_per_kpt
kshift = (k-1)*mo_num_per_kpt
do i = 1, mo_num_per_kpt
if(mo_class(i+kshift) == 'Active')then
n_act_orb_kpts(k) += 1
@ -511,7 +511,7 @@ BEGIN_PROVIDER [ integer, n_virt_orb_kpts, (kpt_num)]
do k=1,kpt_num
n_virt_orb_kpts(k) = 0
kshift = (1-k)*mo_num_per_kpt
kshift = (k-1)*mo_num_per_kpt
do i = 1, mo_num_per_kpt
if(mo_class(i+kshift) == 'Virtual')then
n_virt_orb_kpts(k) += 1
@ -532,7 +532,7 @@ BEGIN_PROVIDER [ integer, n_del_orb_kpts, (kpt_num)]
do k=1,kpt_num
n_del_orb_kpts(k) = 0
kshift = (1-k)*mo_num_per_kpt
kshift = (k-1)*mo_num_per_kpt
do i = 1, mo_num_per_kpt
if(mo_class(i+kshift) == 'Deleted')then
n_del_orb_kpts(k) += 1

101
src/green/EZFIO.cfg Normal file
View File

@ -0,0 +1,101 @@
[n_lanczos_complete]
type: integer
doc: number of lanczos iterations completed
interface: ezfio,provider,ocaml
default: 0
[n_lanczos_iter]
type: integer
doc: number of lanczos iterations
interface: ezfio,provider,ocaml
default: 10
[omega_min]
type: double precision
doc: lower limit of frequency for spectral density calculation
interface: ezfio,provider,ocaml
default: -2.e-1
[omega_max]
type: double precision
doc: upper limit of frequency for spectral density calculation
interface: ezfio,provider,ocaml
default: 1.2e1
[n_omega]
type: integer
doc: number of points for spectral density calculation
interface: ezfio,provider,ocaml
default: 1000
[gf_epsilon]
type: double precision
doc: infinitesimal imaginary frequency term in Green's function
interface: ezfio,provider,ocaml
default: 1.e-2
[n_green_vec]
type: integer
doc: number of holes/particles
interface: ezfio
default: 2
[green_idx]
interface: ezfio
doc: homo/lumo indices
type: integer
size: (green.n_green_vec)
[green_spin]
interface: ezfio
doc: homo/lumo spin
type: integer
size: (green.n_green_vec)
[green_sign]
interface: ezfio
doc: homo/lumo sign
type: double precision
size: (green.n_green_vec)
[alpha_lanczos]
interface: ezfio
doc: lanczos alpha values
type: double precision
size: (green.n_lanczos_iter,green.n_green_vec)
[beta_lanczos]
interface: ezfio
doc: lanczos beta values
type: double precision
size: (green.n_lanczos_iter,green.n_green_vec)
[un_lanczos]
interface: ezfio
doc: saved lanczos u vector
type: double precision
size: (2,determinants.n_det,green.n_green_vec)
[vn_lanczos]
interface: ezfio
doc: saved lanczos v vector
type: double precision
size: (2,determinants.n_det,green.n_green_vec)
[lanczos_eigvals]
interface: ezfio
doc: eigvals of tridiagonal form of H
type: double precision
size: (green.n_lanczos_iter,green.n_green_vec)
[lanczos_debug_print]
interface: ezfio,provider,ocaml
type: logical
doc: printing of lanczos vectors at every step
default: False
[n_lanczos_debug]
interface: ezfio,provider,ocaml
type: integer
doc: number of elements to print
default: 10

1
src/green/NEED Normal file
View File

@ -0,0 +1 @@
davidson fci

6
src/green/README.rst Normal file
View File

@ -0,0 +1,6 @@
=====
dummy
=====
Module necessary to avoid the ``xxx is a root module but does not contain a main file`` message.

51
src/green/green.main.irp.f Normal file
View File

@ -0,0 +1,51 @@
program green
implicit none
BEGIN_DOC
! TODO
END_DOC
read_wf = .True.
touch read_wf
provide n_green_vec
print*,'ref_bitmask_energy = ',ref_bitmask_energy
! call psicoefprinttest
call print_lanczos_eigvals
call print_spec
end
subroutine psicoefprinttest
implicit none
integer :: i
TOUCH psi_coef_complex
print *, 'printing ndet', N_det
end
subroutine print_lanczos_eigvals
implicit none
integer :: i, iunit, j
integer :: getunitandopen
character(5) :: jstr
do j=1,n_green_vec
write(jstr,'(I0.3)') j
iunit = getunitandopen('lanczos_eigval_alpha_beta.out.'//trim(jstr),'w')
print *, 'printing lanczos eigenvalues, alpha, beta to "lanczos_eigval_alpha_beta.out.'//trim(jstr)//'"'
do i=1,n_lanczos_iter
write(iunit,'(I6,3(E25.15))') i, lanczos_eigvals(i,j), alpha_lanczos(i,j), beta_lanczos(i,j)
enddo
close(iunit)
enddo
end
subroutine print_spec
implicit none
integer :: i, iunit, j
integer :: getunitandopen
character(5) :: jstr
do j=1,n_green_vec
write(jstr,'(I0.3)') j
iunit = getunitandopen('omega_A.out.'//trim(jstr),'w')
print *, 'printing frequency, spectral density to "omega_A.out.'//trim(jstr)//'"'
do i=1,n_omega
write(iunit,'(2(E25.15))') omega_list(i), spectral_lanczos(i,j)
enddo
close(iunit)
enddo
end

847
src/green/hu0_hp.irp.f Normal file
View File

@ -0,0 +1,847 @@
! modified from H_S2_u_0_nstates_openmp in Davidson/u0Hu0.irp.f
subroutine h_u_0_hp_openmp(v_0,u_0,N_hp,sze,spin_hp,sign_hp,idx_hp)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
!
! N_hp is number of holes and particles to be applied
! each element of spin_hp is either 1(alpha) or 2(beta)
! each element of sign_hp is either 1(particle) or -1(hole)
! idx_hp contains orbital indices for holes and particles
END_DOC
integer, intent(in) :: N_hp,sze
complex*16, intent(inout) :: v_0(sze,N_hp), u_0(sze,N_hp)
integer :: k
complex*16, allocatable :: u_t(:,:), v_t(:,:)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_hp,N_det),v_t(N_hp,N_det))
do k=1,N_hp
call cdset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
v_t = (0.d0,0.d0)
call cdtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_hp)
call h_u_0_hp_openmp_work(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,1,N_det,0,1)
deallocate(u_t)
call cdtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_hp, N_det)
deallocate(v_t)
do k=1,N_hp
call cdset_order(v_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call cdset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine h_u_0_hp_openmp_work(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_hp,sze,istart,iend,ishift,istep
complex*16, intent(in) :: u_t(N_hp,N_det)
complex*16, intent(out) :: v_t(N_hp,sze)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
PROVIDE ref_bitmask_energy N_int
select case (N_int)
case (1)
call H_u_0_hp_openmp_work_1(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
case (2)
call H_u_0_hp_openmp_work_2(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
case (3)
call H_u_0_hp_openmp_work_3(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
case (4)
call H_u_0_hp_openmp_work_4(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
case default
call H_u_0_hp_openmp_work_N_int(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine h_u_0_hp_openmp_work_$N_int(v_t,u_t,N_hp,sze,spin_hp,sign_hp,idx_hp,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t> and s_t = S^2 |u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_hp,sze,istart,iend,ishift,istep
complex*16, intent(in) :: u_t(N_hp,N_det)
complex*16, intent(out) :: v_t(N_hp,sze)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
complex*16 :: hij
double precision :: hii
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
integer*8 :: k8
logical, allocatable :: exc_is_banned_a1(:),exc_is_banned_b1(:),exc_is_banned_a2(:),exc_is_banned_b2(:)
logical, allocatable :: exc_is_banned_ab1(:),exc_is_banned_ab12(:),allowed_hp(:)
logical :: all_banned_a1,all_banned_b1,all_banned_a2,all_banned_b2
logical :: all_banned_ab12,all_banned_ab1
integer :: ii,na,nb
double precision, allocatable :: hii_hp(:)
complex*16, allocatable :: hij_hp(:)
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson elec_num
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int, &
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_hp, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc, &
!$OMP istart, iend, istep, irp_here, v_t, &
!$OMP spin_hp,sign_hp,idx_hp, &
!$OMP elec_num_tab,nuclear_repulsion, &
!$OMP ishift, idx0, u_t, maxab) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, hii, hij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP exc_is_banned_a1,exc_is_banned_b1,exc_is_banned_ab1, &
!$OMP exc_is_banned_a2,exc_is_banned_b2,exc_is_banned_ab12, &
!$OMP all_banned_a1,all_banned_b1,all_banned_ab1, &
!$OMP all_banned_a2,all_banned_b2,all_banned_ab12, &
!$OMP allowed_hp, &
!$OMP ii, hij_hp, j, hii_hp,na,nb, &
!$OMP n_singles_b, k8)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab), &
exc_is_banned_a1(N_hp), &
exc_is_banned_b1(N_hp), &
exc_is_banned_a2(N_hp), &
exc_is_banned_b2(N_hp), &
exc_is_banned_ab1(N_hp), &
exc_is_banned_ab12(N_hp), &
allowed_hp(N_hp), &
hij_hp(N_hp), &
hii_hp(N_hp))
kcol_prev=-1
all_banned_b1=.False.
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! iterate over dets in psi
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then !if we've moved to a new unique beta determinant
call get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned_b1,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b1)
if (all_banned_b1) then
kcol_prev = kcol
cycle
else ! get all unique beta dets connected to this one by a single excitation
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
kcol_prev = kcol
endif
else
if (all_banned_b1) cycle
endif
! at least some beta allowed
! check alpha
call get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned_a1,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a1)
if (all_banned_a1) cycle
all_banned_ab1=.True.
do ii=1,N_hp
exc_is_banned_ab1(ii)=(exc_is_banned_a1(ii).or.exc_is_banned_b1(ii))
all_banned_ab1 = (all_banned_ab1.and.exc_is_banned_ab1(ii))
enddo
if (all_banned_ab1) cycle
! kcol_prev = kcol ! keep track of old col to see when we've moved to a new one
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b ! loop over other columns in this row
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2)
if (all_banned_b2) cycle
l_a = psi_bilinear_matrix_columns_loc(lcol) ! location of start of this column within psi_bilinear_mat
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a ! loop over rows in this column
lrow = psi_bilinear_matrix_rows(l_a) ! get row (index of unique alpha det)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) ! get alpha det
ASSERT (l_a <= N_det)
idx(j) = l_a ! indices of dets within psi_bilinear_mat
l_a = l_a+1
enddo
j = j-1
! get all alpha dets in this column that are connected to ref alpha by a single exc.
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2)
if (all_banned_a2) cycle
all_banned_ab12 = .True.
do ii=1,N_hp
exc_is_banned_ab12(ii)=((exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii)).or.exc_is_banned_a2(ii))
allowed_hp(ii)=(.not.exc_is_banned_ab12(ii))
all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii))
enddo
if (all_banned_ab12) cycle
call i_h_j_double_alpha_beta_hp(tmp_det,tmp_det2,$N_int,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
do l=1,N_hp
v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1)
if (all_banned_ab1) cycle
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
!call get_all_spin_singles_and_doubles_$N_int( &
! buffer, idx, spindet, i, &
! singles_a, doubles, n_singles_a, n_doubles )
call get_all_spin_singles_and_doubles( &
buffer, idx, spindet, $N_int, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2)
if (all_banned_a2) cycle
all_banned_ab12 = .True.
do ii=1,N_hp
exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_a2(ii))
allowed_hp(ii)=(.not.exc_is_banned_ab12(ii))
all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii))
enddo
if (all_banned_ab12) cycle
call i_h_j_mono_spin_hp(tmp_det,tmp_det2,$N_int,1, hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
do l=1,N_hp
v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
call get_list_hp_banned_single_spin(psi_det_alpha_unique(1,lrow),N_hp,exc_is_banned_a2,spin_hp,sign_hp,idx_hp,1,$N_int,all_banned_a2)
if (all_banned_a2) cycle
all_banned_ab12 = .True.
do ii=1,N_hp
exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_a2(ii))
allowed_hp(ii)=(.not.exc_is_banned_ab12(ii))
all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii))
enddo
if (all_banned_ab12) cycle
call i_h_j_double_spin_hp( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int,1,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
do l=1,N_hp
v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
!! should already be done from top of loop?
!call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1)
!if (all_banned_ab1) cycle
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
!call get_all_spin_singles_and_doubles_$N_int( &
! buffer, idx, spindet, i, &
! singles_b, doubles, n_singles_b, n_doubles )
call get_all_spin_singles_and_doubles( &
buffer, idx, spindet, $N_int, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
call get_list_hp_banned_spin(tmp_det2,N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2)
if (all_banned_b2) cycle
all_banned_ab12 = .True.
do ii=1,N_hp
exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii))
allowed_hp(ii)=(.not.exc_is_banned_ab12(ii))
all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii))
enddo
if (all_banned_ab12) cycle
call i_h_j_mono_spin_hp(tmp_det,tmp_det2,$N_int,2, hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_hp
v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all beta doubles
! ----------------------------------
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
call get_list_hp_banned_single_spin(psi_det_beta_unique(1,lcol),N_hp,exc_is_banned_b2,spin_hp,sign_hp,idx_hp,2,$N_int,all_banned_b2)
if (all_banned_b2) cycle
all_banned_ab12 = .True.
do ii=1,N_hp
exc_is_banned_ab12(ii)=(exc_is_banned_ab1(ii).or.exc_is_banned_b2(ii))
allowed_hp(ii)=(.not.exc_is_banned_ab12(ii))
all_banned_ab12 = (all_banned_ab12.and.exc_is_banned_ab12(ii))
enddo
if (all_banned_ab12) cycle
call i_h_j_double_spin_hp( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int,2,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_hp
v_t(l,k_a) = v_t(l,k_a) + hij_hp(l) * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
call get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned_ab1,spin_hp,sign_hp,idx_hp,$N_int,all_banned_ab1)
if (all_banned_ab1) cycle
double precision, external :: diag_H_mat_elem, diag_S_mat_elem
hii = diag_h_mat_elem(tmp_det,$N_int)
do ii=1,N_hp
if(exc_is_banned_ab1(ii)) then
hii_hp(ii)=0.d0
else
tmp_det2=tmp_det
na=elec_num_tab(spin_hp(ii))
nb=elec_num_tab(iand(spin_hp(ii),1)+1)
hii_hp(ii)=hii
if (sign_hp(ii)>0) then
call ac_operator(idx_hp(ii),spin_hp(ii),tmp_det2,hii_hp(ii),$N_int,na,nb)
else
call a_operator(idx_hp(ii),spin_hp(ii),tmp_det2,hii_hp(ii),$N_int,na,nb)
endif
endif
v_t(ii,k_a) = v_t(ii,k_a) + (nuclear_repulsion + hii_hp(ii)) * u_t(ii,k_a)
enddo
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx, &
exc_is_banned_a1, &
exc_is_banned_b1, &
exc_is_banned_a2, &
exc_is_banned_b2, &
exc_is_banned_ab1, &
exc_is_banned_ab12, &
allowed_hp, &
hij_hp, hii_hp )
!$OMP END PARALLEL
deallocate(idx0)
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE
subroutine i_h_j_double_spin_hp(key_i,key_j,Nint,ispin,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
use bitmasks
implicit none
BEGIN_DOC
! todo: maybe make new get_double_excitation_spin?
! the 4 index ordering is already done in there, so we could avoid duplicating that work
! Returns <i|H|j> where i and j are determinants differing by a same-spin double excitation
END_DOC
integer, intent(in) :: Nint,ispin,N_hp
integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint)
complex*16, intent(out) :: hij_hp(N_hp)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
logical, intent(in) :: allowed_hp(N_hp)
complex*16 :: hij0
double precision :: phase_hp(N_hp)
integer :: exc(0:2,2)
double precision :: phase
complex*16, external :: mo_two_e_integral_complex
integer :: i1,i2,i3,i4,j2,j3,ii
PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map
call get_double_excitation_spin(key_i,key_j,exc,phase,Nint)
hij0 = phase*(mo_two_e_integral_complex( &
exc(1,1), &
exc(2,1), &
exc(1,2), &
exc(2,2)) - &
mo_two_e_integral_complex( &
exc(1,1), &
exc(2,1), &
exc(2,2), &
exc(1,2)) )
ASSERT (exc(1,1) < exc(2,1))
ASSERT (exc(1,2) < exc(2,2))
i1=min(exc(1,1),exc(1,2))
j2=max(exc(1,1),exc(1,2))
j3=min(exc(2,1),exc(2,2))
i4=max(exc(2,1),exc(2,2))
i2=min(j2,j3)
i3=max(j2,j3)
do ii=1,N_hp
if (allowed_hp(ii)) then
if (ispin.eq.spin_hp(ii)) then
if ((idx_hp(ii).lt.i1).or.(idx_hp(ii).gt.i4)) then
phase_hp(ii)=1.d0
else if ((idx_hp(ii).lt.i2).or.(idx_hp(ii).gt.i3)) then
phase_hp(ii)=-1.d0
else
phase_hp(ii)=1.d0
endif
else
phase_hp(ii)=1.d0
endif
else
phase_hp(ii)=0.d0
endif
hij_hp(ii) = hij0 * phase_hp(ii)
enddo
end
subroutine i_h_j_mono_spin_hp(key_i,key_j,Nint,spin,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
use bitmasks
implicit none
BEGIN_DOC
! todo: change this to use normal version of get_mono_excitation_from_fock
! all info needed is in phase and hij, h/p part can happen after getting hij the normal way
! Returns <i|H|j> where i and j are determinants differing by a single excitation
END_DOC
integer, intent(in) :: Nint, spin, N_hp
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
complex*16, intent(out) :: hij_hp(N_hp)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
logical, intent(in) :: allowed_hp(N_hp)
!double precision :: phase_hp(N_hp)
complex*16 :: hij0
integer :: exc(0:2,2)
double precision :: phase
PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map
call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint)
call get_single_excitation_from_fock_hp(key_i,key_j,exc(1,1),exc(1,2),spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
end
subroutine get_single_excitation_from_fock_hp(det_1,det_2,h,p,spin,phase,N_hp,hij_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
use bitmasks
implicit none
integer,intent(in) :: h,p,spin,N_hp
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2)
complex*16, intent(out) :: hij_hp(N_hp)
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
logical, intent(in) :: allowed_hp(N_hp)
double precision :: phase_hp(N_hp)
complex*16 :: hij0
integer :: low,high
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i,ii
do i = 1, N_int
differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
partcl(i,1) = iand(differences(i,1),det_1(i,1))
partcl(i,2) = iand(differences(i,2),det_1(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
hij0 = fock_op_cshell_ref_bitmask_cplx(h,p)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)
hij0 -= big_array_coulomb_integrals_complex(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
hij0 -= big_array_coulomb_integrals_complex(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
hij0 += big_array_exchange_integrals_complex(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
hij0 += big_array_coulomb_integrals_complex(i,h,p)!get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
hij0 += big_array_coulomb_integrals_complex(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
hij0 -= big_array_exchange_integrals_complex(i,h,p)!get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map)
enddo
low=min(h,p)
high=max(h,p)
!! do ii=1,N_hp
!! if (.not.allowed_hp(ii)) then
!! phase_hp(ii) = 0.d0
!! cycle
!! else if (spin_hp(ii).ne.spin) then
!! phase_hp(ii) = 1.d0
!! else
!! if ((low.lt.idx_hp(ii)).and.(high.gt.idx_hp(ii))) then
!! phase_hp(ii) = -1.d0
!! else
!! phase_hp(ii) = 1.d0
!! endif
!! endif
!! enddo
!!
!! do ii=1,N_hp
!! if (allowed_hp(ii)) then
!! hij_hp(ii) = hij + sign_hp(ii) * big_array_coulomb_integrals(idx_hp(ii),h,p)
!! if (spin.eq.spin_hp(ii)) then
!! hij_hp(ii) = hij_hp(ii) - sign_hp(ii) * big_array_exchange_integrals(idx_hp(ii),h,p)
!! endif
!! else
!! hij_hp(ii) = 0.d0
!! endif
!! enddo
!!
!! do ii=1,N_hp
!! hij_hp(ii) = hij_hp(ii) * phase_hp(ii) * phase
!! enddo
do ii=1,N_hp
if (.not.allowed_hp(ii)) then
phase_hp(ii) = 0.d0
hij_hp(ii) = 0.d0
cycle
else if (spin.eq.spin_hp(ii)) then
hij_hp(ii) = hij0 + sign_hp(ii) *(big_array_coulomb_integrals_complex(idx_hp(ii),h,p) - big_array_exchange_integrals_complex(idx_hp(ii),h,p))
if ((low.lt.idx_hp(ii)).and.(high.gt.idx_hp(ii))) then
phase_hp(ii) = -1.d0
else
phase_hp(ii) = 1.d0
endif
else
phase_hp(ii) = 1.d0
hij_hp(ii) = hij0 + sign_hp(ii) * big_array_coulomb_integrals_complex(idx_hp(ii),h,p)
endif
hij_hp(ii) = hij_hp(ii) * phase * phase_hp(ii)
enddo
end
subroutine i_H_j_double_alpha_beta_hp(key_i,key_j,Nint,hij_hp,N_hp,spin_hp,sign_hp,idx_hp,allowed_hp)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants differing by an opposite-spin double excitation
END_DOC
integer, intent(in) :: Nint,N_hp
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
complex*16, intent(out) :: hij_hp(N_hp)
complex*16 :: hij0
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
logical, intent(in) :: allowed_hp(N_hp)
double precision :: phase_hp(N_hp)
integer :: i
integer :: lowhigh(2,2)
integer :: exc(0:2,2,2)
double precision :: phase, phase2
complex*16, external :: mo_two_e_integral_complex
PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map
call get_single_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint)
call get_single_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint)
phase = phase*phase2
if (exc(1,1,1) == exc(1,2,2)) then
hij0 = big_array_exchange_integrals_complex(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) == exc(1,1,2)) then
hij0 = big_array_exchange_integrals_complex(exc(1,2,1),exc(1,1,1),exc(1,2,2))
else
hij0 = mo_two_e_integral_complex( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2))
endif
!todo: clean this up
! if new particle/hole is between p/h of single exc of same spin, then parity changes, otherwise stays the same
! value of Hij for double excitation is unchanged (new p/h is not one of the indices involved in the excitation)
lowhigh(1,1)=min(exc(1,1,1),exc(1,2,1))
lowhigh(2,1)=max(exc(1,1,1),exc(1,2,1))
lowhigh(1,2)=min(exc(1,1,2),exc(1,2,2))
lowhigh(2,2)=max(exc(1,1,2),exc(1,2,2))
do i=1,N_hp
if (.not.allowed_hp(i)) then
phase_hp(i)=0.d0
else if ((idx_hp(i).gt.lowhigh(1,spin_hp(i))).and.(idx_hp(i).lt.lowhigh(2,spin_hp(i)))) then
phase_hp(i)=-1.d0
else
phase_hp(i)=1.d0
endif
hij_hp(i)=hij0*phase*phase_hp(i)
enddo
end

405
src/green/hu0_lanczos.irp.f Normal file
View File

@ -0,0 +1,405 @@
! modified from H_S2_u_0_nstates_openmp in Davidson/u0Hu0.irp.f
subroutine H_u_0_openmp(v_0,u_0,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer :: N_st=1
integer, intent(in) :: sze
complex*16, intent(inout) :: v_0(sze), u_0(sze)
integer :: k
call cdset_order(u_0(1),psi_bilinear_matrix_order,N_det)
v_0 = (0.d0,0.d0)
call h_u_0_openmp_work(v_0,u_0,sze,1,N_det,0,1)
call cdset_order(v_0(1),psi_bilinear_matrix_order_reverse,N_det)
call cdset_order(u_0(1),psi_bilinear_matrix_order_reverse,N_det)
end
subroutine H_u_0_openmp_work(v_t,u_t,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer :: N_st=1
integer, intent(in) :: sze,istart,iend,ishift,istep
complex*16, intent(in) :: u_t(N_det)
complex*16, intent(out) :: v_t(sze)
PROVIDE ref_bitmask_energy N_int
select case (N_int)
case (1)
call H_u_0_openmp_work_1(v_t,u_t,sze,istart,iend,ishift,istep)
case (2)
call H_u_0_openmp_work_2(v_t,u_t,sze,istart,iend,ishift,istep)
case (3)
call H_u_0_openmp_work_3(v_t,u_t,sze,istart,iend,ishift,istep)
case (4)
call H_u_0_openmp_work_4(v_t,u_t,sze,istart,iend,ishift,istep)
case default
call H_u_0_openmp_work_N_int(v_t,u_t,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine H_u_0_openmp_work_$N_int(v_t,u_t,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer :: N_st=1
integer, intent(in) :: sze,istart,iend,ishift,istep
complex*16, intent(in) :: u_t(N_det)
complex*16, intent(out) :: v_t(sze)
complex*16 :: hij
double precision :: hii
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
integer*8 :: k8
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int, &
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc, &
!$OMP istart, iend, istep, irp_here, v_t, &
!$OMP ishift, idx0, u_t, maxab) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, hii, hij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, k8)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_h_j_double_alpha_beta_complex(tmp_det,tmp_det2,$N_int,hij)
v_t(k_a) = v_t(k_a) + hij * u_t(l_a)
enddo
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
!call get_all_spin_singles_and_doubles_$N_int( &
! buffer, idx, spindet, i, &
! singles_a, doubles, n_singles_a, n_doubles )
call get_all_spin_singles_and_doubles( &
buffer, idx, spindet, $N_int, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_h_j_single_spin_complex( tmp_det, tmp_det2, $N_int, 1, hij)
v_t(k_a) = v_t(k_a) + hij * u_t(l_a)
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
call i_h_j_double_spin_complex( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
v_t(k_a) = v_t(k_a) + hij * u_t(l_a)
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
!call get_all_spin_singles_and_doubles_$N_int( &
! buffer, idx, spindet, i, &
! singles_b, doubles, n_singles_b, n_doubles )
call get_all_spin_singles_and_doubles( &
buffer, idx, spindet, $N_int, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
call i_h_j_single_spin_complex( tmp_det, tmp_det2, $N_int, 2, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
v_t(k_a) = v_t(k_a) + hij * u_t(l_a)
enddo
! Compute Hij for all beta doubles
! ----------------------------------
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
call i_h_j_double_spin_complex( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
v_t(k_a) = v_t(k_a) + hij * u_t(l_a)
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
!double precision, external :: diag_H_mat_elem, diag_S_mat_elem
double precision, external :: diag_H_mat_elem
hii = diag_H_mat_elem(tmp_det,$N_int)
v_t(k_a) = v_t(k_a) + hii * u_t(k_a)
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE

883
src/green/lanczos.irp.f Normal file
View File

@ -0,0 +1,883 @@
BEGIN_PROVIDER [ integer, n_green_vec ]
implicit none
BEGIN_DOC
! number of particles/holes to use for spectral density calc.
! just set to 2 for now (homo and lumo)
END_DOC
n_green_vec = 2
END_PROVIDER
BEGIN_PROVIDER [ integer, green_idx, (n_green_vec) ]
&BEGIN_PROVIDER [ integer, green_idx_int, (n_green_vec) ]
&BEGIN_PROVIDER [ integer, green_idx_bit, (n_green_vec) ]
&BEGIN_PROVIDER [ integer, green_spin, (n_green_vec) ]
&BEGIN_PROVIDER [ double precision, green_sign, (n_green_vec) ]
implicit none
BEGIN_DOC
! description of particles/holes to be used in spectral density calculation
! green_idx: orbital index of particle/hole
! green_idx_{int,bit}: location of idx within determinant bitstring
! green_spin: 1(alpha) or 2(beta)
! green_sign: 1(particle) or -1(hole)
END_DOC
integer :: s1,s2,i1,i2
integer :: i
integer :: idx_homo_lumo(2), spin_homo_lumo(2)
logical :: has_idx,has_spin,has_sign,has_lanc
integer :: nlanc
! needs psi_det, mo_num, N_int, mo_bielec_integral_jj, mo_mono_elec_integral_diag
call ezfio_has_green_green_idx(has_idx)
call ezfio_has_green_green_spin(has_spin)
call ezfio_has_green_green_sign(has_sign)
! call ezfio_has_green_n_lanczos_complete(has_lanc)
call ezfio_get_green_n_lanczos_complete(nlanc)
if (has_idx.and.has_spin.and.has_sign) then
print*,'reading idx,spin,sign'
call ezfio_get_green_green_idx(green_idx)
call ezfio_get_green_green_spin(green_spin)
call ezfio_get_green_green_sign(green_sign)
else if (nlanc.gt.0) then
stop 'problem with lanczos restart; need idx, spin, sign'
else
print*,'new lanczos calculation, finding homo/lumo'
call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_num,idx_homo_lumo,spin_homo_lumo)
! homo
green_idx(1)=idx_homo_lumo(1)
green_spin(1)=spin_homo_lumo(1)
green_sign(1)=-1.d0
! lumo
green_idx(2)=idx_homo_lumo(2)
green_spin(2)=spin_homo_lumo(2)
green_sign(2)=1.d0
call ezfio_set_green_green_idx(green_idx)
call ezfio_set_green_green_spin(green_spin)
call ezfio_set_green_green_sign(green_sign)
endif
! if (nlanc.gt.0) then
! ! call ezfio_get_green_n_lanczos_complete(nlanc)
! print*,'restarting from previous lanczos',nlanc
! if (has_idx.and.has_spin.and.has_sign) then
! print*,'reading idx,spin,sign'
! call ezfio_get_green_green_idx(green_idx)
! call ezfio_get_green_green_spin(green_spin)
! call ezfio_get_green_green_sign(green_sign)
! else
! stop 'problem with lanczos restart; need idx, spin, sign'
! endif
! else
! print*,'new lanczos calculation, finding homo/lumo'
! call get_homo_lumo(psi_det(1:N_int,1:2,1),N_int,mo_num,idx_homo_lumo,spin_homo_lumo)
!
! ! homo
! green_idx(1)=idx_homo_lumo(1)
! green_spin(1)=spin_homo_lumo(1)
! green_sign(1)=-1.d0
!
! ! lumo
! green_idx(2)=idx_homo_lumo(2)
! green_spin(2)=spin_homo_lumo(2)
! green_sign(2)=1.d0
!
! call ezfio_set_green_green_idx(green_idx)
! call ezfio_set_green_green_spin(green_spin)
! call ezfio_set_green_green_sign(green_sign)
! endif
do i=1,n_green_vec
call get_orb_int_bit(green_idx(i),green_idx_int(i),green_idx_bit(i))
print*,i,green_idx(i),green_idx_int(i),green_idx_bit(i),green_spin(i),green_sign(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, green_det_phase, (N_det,n_green_vec) ]
implicit none
BEGIN_DOC
! for each det in psi, compute phase for each particle/hole excitation
! each element should be +/-1 or 0
END_DOC
integer :: i
double precision :: phase_tmp(n_green_vec)
PROVIDE psi_det green_idx
do i=1,N_det
call get_phase_hp(green_idx_int,green_idx_bit,green_spin,green_sign,psi_det(1,1,i),phase_tmp,N_int,n_green_vec)
green_det_phase(i,1:n_green_vec) = phase_tmp(1:n_green_vec)
enddo
END_PROVIDER
BEGIN_PROVIDER [ complex*16, u1_lanczos, (N_det,n_green_vec) ]
implicit none
BEGIN_DOC
! initial lanczos vectors
! must be normalized
END_DOC
integer :: i,j
do j=1,n_green_vec
do i=1,N_det
u1_lanczos(i,j)=green_det_phase(i,j)*psi_coef(i,1)
enddo
call normalize_complex(u1_lanczos(:,j),N_det)
enddo
END_PROVIDER
! BEGIN_PROVIDER [ double precision, alpha_lanczos, (n_green_vec,n_lanczos_iter) ]
!&BEGIN_PROVIDER [ double precision, beta_lanczos, (n_green_vec,n_lanczos_iter) ]
BEGIN_PROVIDER [ double precision, alpha_lanczos, (n_lanczos_iter,n_green_vec) ]
&BEGIN_PROVIDER [ double precision, beta_lanczos, (n_lanczos_iter,n_green_vec) ]
&BEGIN_PROVIDER [ complex*16, un_lanczos, (N_det,n_green_vec) ]
&BEGIN_PROVIDER [ complex*16, vn_lanczos, (N_det,n_green_vec) ]
&BEGIN_PROVIDER [ double precision, lanczos_eigvals, (n_lanczos_iter,n_green_vec) ]
implicit none
BEGIN_DOC
! for each particle/hole:
! provide alpha and beta for tridiagonal form of H
! un, vn lanczos vectors from latest iteration
! lanczos_eigvals: eigenvalues of tridiagonal form of H
END_DOC
PROVIDE lanczos_debug_print n_lanczos_debug
complex*16, allocatable :: work(:,:)
! double precision :: alpha_tmp,beta_tmp
double precision, allocatable :: alpha_tmp(:),beta_tmp(:)
double precision, allocatable :: alpha_tmp_vec(:,:), beta_tmp_vec(:,:)
integer :: i,j
integer :: n_lanc_new_tmp, n_lanc_old_tmp
call ezfio_get_green_n_lanczos_iter(n_lanc_new_tmp)
call ezfio_get_green_n_lanczos_complete(n_lanc_old_tmp)
if ((n_lanczos_complete).gt.0) then
! allocate(alpha_tmp_vec(n_green_vec,n_lanczos_complete),beta_tmp_vec(n_green_vec,n_lanczos_complete))
allocate(alpha_tmp_vec(n_lanczos_complete,n_green_vec),beta_tmp_vec(n_lanczos_complete,n_green_vec))
logical :: has_un_lanczos, has_vn_lanczos
call ezfio_has_green_un_lanczos(has_un_lanczos)
call ezfio_has_green_vn_lanczos(has_vn_lanczos)
if (has_un_lanczos.and.has_vn_lanczos) then
call ezfio_get_green_un_lanczos(un_lanczos)
call ezfio_get_green_vn_lanczos(vn_lanczos)
! if (lanczos_debug_print) then
! print*,'uu,vv read from disk'
! do i=1,n_lanczos_debug
! write(6,'(4(E25.15))')un_lanczos(i),vn_lanczos(i)
! enddo
! endif
else
print*,'problem reading lanczos vectors for restart'
stop
endif
logical :: has_alpha_lanczos, has_beta_lanczos
call ezfio_has_green_alpha_lanczos(has_alpha_lanczos)
call ezfio_has_green_beta_lanczos(has_beta_lanczos)
if (has_alpha_lanczos.and.has_beta_lanczos) then
call ezfio_set_green_n_lanczos_iter(n_lanc_old_tmp)
call ezfio_get_green_alpha_lanczos(alpha_tmp_vec)
call ezfio_get_green_beta_lanczos(beta_tmp_vec)
call ezfio_set_green_n_lanczos_iter(n_lanc_new_tmp)
do j=1,n_green_vec
do i=1,n_lanczos_complete
alpha_lanczos(i,j)=alpha_tmp_vec(i,j)
beta_lanczos(i,j)=beta_tmp_vec(i,j)
enddo
enddo
else
print*,'problem reading lanczos alpha, beta for restart'
stop
endif
deallocate(alpha_tmp_vec,beta_tmp_vec)
else
call write_time(6)
print*,'no saved lanczos vectors. starting lanczos'
PROVIDE u1_lanczos
un_lanczos=u1_lanczos
allocate(work(N_det,n_green_vec),alpha_tmp(n_green_vec),beta_tmp(n_green_vec))
call lanczos_h_init_hp(un_lanczos,vn_lanczos,work,N_det,alpha_tmp,beta_tmp,&
n_green_vec,green_spin,green_sign,green_idx)
do i=1,n_green_vec
alpha_lanczos(1,i)=alpha_tmp(i)
beta_lanczos(1,i)=beta_tmp(i)
enddo
n_lanczos_complete=1
deallocate(work,alpha_tmp,beta_tmp)
endif
allocate(work(N_det,n_green_vec),alpha_tmp(n_green_vec),beta_tmp(n_green_vec))
do i=n_lanczos_complete+1,n_lanczos_iter
call write_time(6)
print*,'starting lanczos iteration',i
call lanczos_h_step_hp(un_lanczos,vn_lanczos,work,N_det,alpha_tmp,beta_tmp,&
n_green_vec,green_spin,green_sign,green_idx)
do j=1,n_green_vec
alpha_lanczos(i,j)=alpha_tmp(j)
beta_lanczos(i,j)=beta_tmp(j)
enddo
n_lanczos_complete=n_lanczos_complete+1
enddo
deallocate(work,alpha_tmp,beta_tmp)
call ezfio_set_green_alpha_lanczos(alpha_lanczos)
call ezfio_set_green_beta_lanczos(beta_lanczos)
call ezfio_set_green_un_lanczos(un_lanczos)
call ezfio_set_green_vn_lanczos(vn_lanczos)
call ezfio_set_green_n_lanczos_complete(n_lanczos_complete)
call diag_lanczos_vals_hp(alpha_lanczos, beta_lanczos, n_lanczos_complete, lanczos_eigvals,&
n_lanczos_iter,n_green_vec)
call ezfio_set_green_lanczos_eigvals(lanczos_eigvals)
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_omega ]
implicit none
BEGIN_DOC
! step size between frequency points for spectral density calculation
! calculated from min, max, and number of steps
END_DOC
delta_omega=(omega_max-omega_min)/n_omega
END_PROVIDER
BEGIN_PROVIDER [ double precision, omega_list, (n_omega) ]
implicit none
BEGIN_DOC
! list of frequencies at which to compute spectral density
END_DOC
integer :: i
double precision :: omega_i
PROVIDE delta_omega
do i=1,n_omega
omega_list(i) = omega_min + (i-1)*delta_omega
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, spectral_lanczos, (n_omega,n_green_vec) ]
implicit none
BEGIN_DOC
! spectral density A(omega) calculated from lanczos alpha/beta
! calculated for n_omega points between omega_min and omega_max
END_DOC
integer :: i,j
double precision :: omega_i
complex*16 :: z_i
!double precision :: spec_lanc_rev
double precision :: spec_lanc_rev_sign
logical :: has_ci_energy
double precision :: ref_energy_0
PROVIDE delta_omega alpha_lanczos beta_lanczos omega_list
call ezfio_has_fci_energy(has_ci_energy)
if (has_ci_energy) then
call ezfio_get_fci_energy(ref_energy_0)
else
print*,'no reference energy from full_ci_zmq, exiting'
stop
endif
do i=1,n_omega
omega_i = omega_list(i)
z_i = dcmplx(omega_i,gf_epsilon)
do j=1,n_green_vec
! spectral_lanczos(i,j) = spec_lanc_rev(n_lanczos_iter,alpha_lanczos(:,j),beta_lanczos(:,j),z_i)
spectral_lanczos(i,j) = spec_lanc_rev_sign(n_lanczos_iter, &
alpha_lanczos(:,j), &
beta_lanczos(:,j), &
z_i - green_sign(j)*ref_energy_0, &
green_sign(j))
enddo
enddo
END_PROVIDER
double precision function spec_lanc(n_lanc_iter,alpha,beta,z)
include 'constants.include.F'
implicit none
BEGIN_DOC
! input:
! alpha, beta: from tridiagonal form of H (obtain via lanczos)
! beta and alpha same size (beta(1) is not used)
! n_lanc_iter: size of alpha, beta
! z: omega + i*epsilon
! omega is frequency for which spectral density is to be computed
! epsilon is magnitude of infinitesimal imaginary term
! output:
! spec_lanc: spectral density A(omega)
!
! uses inv_pi=(1.d0/pi) from constants
END_DOC
integer, intent(in) :: n_lanc_iter
double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter)
complex*16, intent(in) :: z
complex*16 bigAj2,bigAj1,bigAj0
complex*16 bigBj2,bigBj1,bigBj0
integer :: j
! init for j=1
! bigAj2 is A(j-2)
! bigAj1 is A(j-1)
! etc.
bigAj2=1.d0 ! A(-1)
bigAj1=0.d0 ! A(0)
bigAj0=1.d0 ! A(1)
bigBj2=0.d0 ! B(-1)
bigBj1=1.d0 ! B(0)
bigBj0=z-alpha(1) ! B(1)
do j=2,n_lanc_iter
bigAj2=bigAj1
bigAj1=bigAj0
bigAj0=(z-alpha(j))*bigAj1 - beta(j)**2*bigAj2
bigBj2=bigBj1
bigBj1=bigBj0
bigBj0=(z-alpha(j))*bigBj1 - beta(j)**2*bigBj2
enddo
spec_lanc=-imag(bigAj0/bigBj0)*inv_pi
end
double precision function spec_lanc_rev(n_lanc_iter,alpha,beta,z)
include 'constants.include.F'
implicit none
BEGIN_DOC
! reverse iteration is more numerically stable
! input:
! alpha, beta: from tridiagonal form of H (obtain via lanczos)
! beta and alpha same size (beta(1) is not used)
! n_lanc_iter: size of alpha, beta
! z: omega + i*epsilon
! omega is frequency for which spectral density is to be computed
! epsilon is magnitude of infinitesimal imaginary term
! output:
! spec_lanc: spectral density A(omega)
!
! uses inv_pi=(1.d0/pi) from constants
END_DOC
integer, intent(in) :: n_lanc_iter
double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter)
complex*16, intent(in) :: z
complex*16 :: tmp
integer :: j
tmp=(0.d0,0.d0)
do j=n_lanc_iter,2,-1
tmp=-beta(j)**2/(z-alpha(j)+tmp)
enddo
tmp=1.d0/(z-alpha(1)+tmp)
spec_lanc_rev=-imag(tmp)*inv_pi
end
double precision function spec_lanc_rev_sign(n_lanc_iter,alpha,beta,z,g_sign)
include 'constants.include.F'
implicit none
BEGIN_DOC
! reverse iteration is more numerically stable
! input:
! alpha, beta: from tridiagonal form of H (obtain via lanczos)
! beta and alpha same size (beta(1) is not used)
! n_lanc_iter: size of alpha, beta
! z: omega + i*epsilon
! omega is frequency for which spectral density is to be computed
! epsilon is magnitude of infinitesimal imaginary term
! output:
! spec_lanc: spectral density A(omega)
!
! uses inv_pi=(1.d0/pi) from constants
END_DOC
integer, intent(in) :: n_lanc_iter
double precision, intent(in) :: alpha(n_lanc_iter), beta(n_lanc_iter)
complex*16, intent(in) :: z
double precision, intent(in) :: g_sign
complex*16 :: tmp
integer :: j
tmp=(0.d0,0.d0)
do j=n_lanc_iter,2,-1
tmp=-beta(j)**2/(z+g_sign*alpha(j)+tmp)
enddo
tmp=1.d0/(z+g_sign*alpha(1)+tmp)
spec_lanc_rev_sign=-imag(tmp)*inv_pi
end
subroutine lanczos_h_init_hp(uu,vv,work,sze,alpha_i,beta_i,ng,spin_hp,sign_hp,idx_hp)
implicit none
integer, intent(in) :: sze,ng
complex*16, intent(in) :: uu(sze,ng)
complex*16, intent(out) :: vv(sze,ng)
complex*16 :: work(sze,ng)
double precision, intent(out) :: alpha_i(ng), beta_i(ng)
integer, intent(in) :: spin_hp(ng), idx_hp(ng)
double precision, intent(in) :: sign_hp(ng)
double precision, external :: dznrm2
complex*16, external :: u_dot_v_complex
integer :: i,j
BEGIN_DOC
! initial step for lanczos tridiagonalization of H for multiple holes/particles
! uu is array of initial vectors u1 (creation/annihilation operator applied to psi)
! output vv is array of lanczos v1 (one for each hole/particle)
END_DOC
print *,'starting lanczos'
print *,'sze = ',sze
! |uu> is |u(1)>
! |w(1)> = H|u(1)>
! |work> is now |w(1)>
call compute_hu_hp(uu,work,ng,sze,spin_hp,sign_hp,idx_hp)
! alpha(n+1) = <u(n+1)|w(n+1)>
do i=1,ng
alpha_i(i)=real(u_dot_v_complex(uu(1:sze,i),work(1:sze,i),sze))
enddo
do j=1,ng
do i=1,sze
vv(i,j)=work(i,j)-alpha_i(j)*uu(i,j)
! write(6,'(7(E25.15))')uu(i,j),vv(i,j),work(i,j),alpha_i(j)
enddo
enddo
beta_i=0.d0
! |vv> is |v(1)>
! |uu> is |u(1)>
end
subroutine lanczos_h_step_hp(uu,vv,work,sze,alpha_i,beta_i,ng,spin_hp,sign_hp,idx_hp)
implicit none
integer, intent(in) :: sze,ng
complex*16, intent(inout) :: uu(sze,ng),vv(sze,ng)
complex*16, intent(out) :: work(sze,ng)
double precision, intent(out) :: alpha_i(ng), beta_i(ng)
integer, intent(in) :: spin_hp(ng), idx_hp(ng)
double precision, intent(in) :: sign_hp(ng)
double precision, external :: dznrm2
complex*16, external :: u_dot_v_complex
integer :: i,j
complex*16 :: tmp_c16
BEGIN_DOC
! lanczos tridiagonalization of H
! n_lanc_iter is number of lanczos iterations
! u1 is initial lanczos vector
! u1 should be normalized
END_DOC
! |vv> is |v(n)>
! |uu> is |u(n)>
! compute beta(n+1)
do j=1,ng
beta_i(j)=dznrm2(sze,vv(:,j),1)
! |vv> is now |u(n+1)>
call zdscal(sze,(1.d0/beta_i(j)),vv(:,j),1)
enddo
! |w(n+1)> = H|u(n+1)>
! |work> is now |w(n+1)>
call compute_hu_hp(vv,work,ng,sze,spin_hp,sign_hp,idx_hp)
! alpha(n+1) = <u(n+1)|w(n+1)>
do i=1,ng
alpha_i(i)=real(u_dot_v_complex(vv(1:sze,i),work(1:sze,i),sze))
enddo
do j=1,ng
do i=1,sze
tmp_c16=work(i,j)-alpha_i(j)*vv(i,j)-beta_i(j)*uu(i,j)
uu(i,j)=vv(i,j)
vv(i,j)=tmp_c16
enddo
enddo
! |vv> is |v(n+1)>
! |uu> is |u(n+1)>
end
subroutine lanczos_h_init(uu,vv,work,sze,alpha_i,beta_i)
implicit none
integer, intent(in) :: sze
complex*16, intent(inout) :: uu(sze)
complex*16, intent(out) :: vv(sze)
complex*16 :: work(sze)
double precision, intent(out) :: alpha_i, beta_i
double precision, external :: dznrm2
complex*16, external :: u_dot_v_complex
integer :: i
BEGIN_DOC
! lanczos tridiagonalization of H
! n_lanc_iter is number of lanczos iterations
! u1 is initial lanczos vector
! u1 should be normalized
END_DOC
print *,'starting lanczos'
print *,'sze = ',sze
! exit if u1 is not normalized
! beta_norm = dznrm2(h_size,u1,1)
! if (dabs(beta_norm-1.d0) .gt. 1.d-6) then
! print *, 'Error: initial Lanczos vector is not normalized'
! stop -1
! endif
! |uu> is |u(1)>
! |w(1)> = H|u(1)>
! |work> is now |w(1)>
call compute_hu(uu,work,sze)
! alpha(n+1) = <u(n+1)|w(n+1)>
alpha_i=real(u_dot_v_complex(uu,work,sze))
do i=1,sze
vv(i)=work(i)-alpha_i*uu(i)
enddo
beta_i=0.d0
if (lanczos_debug_print) then
print*,'init uu,vv,work'
do i=1,n_lanczos_debug
write(6,'(6(E25.15))')uu(i),vv(i),work(i)
enddo
endif
! |vv> is |v(1)>
! |uu> is |u(1)>
end
subroutine lanczos_h_step(uu,vv,work,sze,alpha_i,beta_i)
implicit none
integer, intent(in) :: sze
complex*16, intent(inout) :: uu(sze),vv(sze)
complex*16, intent(out) :: work(sze)
double precision, intent(out) :: alpha_i, beta_i
double precision, external :: dznrm2
complex*16, external :: u_dot_v_complex
integer :: i
complex*16 :: tmp_c16
BEGIN_DOC
! lanczos tridiagonalization of H
! n_lanc_iter is number of lanczos iterations
! u1 is initial lanczos vector
! u1 should be normalized
END_DOC
! exit if u1 is not normalized
! beta_norm = dznrm2(h_size,u1,1)
! if (dabs(beta_norm-1.d0) .gt. 1.d-6) then
! print *, 'Error: initial Lanczos vector is not normalized'
! stop -1
! endif
! |vv> is |v(n)>
! |uu> is |u(n)>
! compute beta(n+1)
beta_i=dznrm2(sze,vv,1)
if (lanczos_debug_print) then
print*,'uu,vv in'
do i=1,n_lanczos_debug
write(6,'(4(E25.15))')uu(i),vv(i)
enddo
endif
! |vv> is now |u(n+1)>
call zdscal(sze,(1.d0/beta_i),vv,1)
! |w(n+1)> = H|u(n+1)>
! |work> is now |w(n+1)>
call compute_hu(vv,work,sze)
if (lanczos_debug_print) then
print*,'vv,work'
do i=1,n_lanczos_debug
write(6,'(4(E25.15))')vv(i),work(i)
enddo
endif
! alpha(n+1) = <u(n+1)|w(n+1)>
alpha_i=real(u_dot_v_complex(vv,work,sze))
do i=1,sze
tmp_c16=work(i)-alpha_i*vv(i)-beta_i*uu(i)
uu(i)=vv(i)
vv(i)=tmp_c16
enddo
! |vv> is |v(n+1)>
! |uu> is |u(n+1)>
end
subroutine lanczos_h(n_lanc_iter,alpha,beta,u1)
implicit none
integer, intent(in) :: n_lanc_iter
double precision, intent(out) :: alpha(n_lanc_iter), beta(n_lanc_iter)
complex*16, intent(in) :: u1(N_det)
integer :: h_size
double precision :: beta_norm, beta_norm_inv
complex*16, allocatable :: vec1(:), vec2(:), vec3(:)
complex*16 :: vec_tmp
double precision, external :: dznrm2
complex*16, external :: u_dot_v_complex
integer :: i,j,l
h_size=N_det
BEGIN_DOC
! lanczos tridiagonalization of H
! n_lanc_iter is number of lanczos iterations
! u1 is initial lanczos vector
! u1 should be normalized
END_DOC
print *,'starting lanczos'
print *,'h_size = ',h_size
! print *,'initial vector:'
! do i=1,h_size
! print *,u1(i)
! enddo
! exit if u1 is not normalized
beta_norm = dznrm2(h_size,u1,1)
if (dabs(beta_norm-1.d0) .gt. 1.d-6) then
print *, 'Error: initial Lanczos vector is not normalized'
stop -1
endif
allocate(vec1(h_size), &
vec2(h_size), &
vec3(h_size))
do i=1,h_size
vec1(i)=u1(i)
enddo
! |w1> = H|u1>
! |vec2> = H|vec1>
call compute_hu(vec1,vec2,h_size)!! TODO: not implemented
! alpha(1) = <u1|H|u1> = <u1|w1>
! = <vec1|vec2>
alpha(1)=real(u_dot_v_complex(vec1,vec2,h_size))
! |v1> = |w1> - alpha(1)*|u1>
! |vec3> = |vec2> - alpha(1)*|vec1>
do i=1,h_size
vec3(i)=vec2(i)-alpha(1)*vec1(i)
enddo
do j=2,n_lanc_iter
call write_time(6)
print *,'starting lanczos iteration:',j
!! vec1 is |u(j-1)>
!! vec3 is |v(j-1)>
! beta(j) = sqrt(<v(j-1)|v(j-1)>)
beta_norm=dznrm2(h_size,vec3,1)
! TODO: check for beta=0?
beta_norm_inv=1.d0/beta_norm
! normalize |v(j-1)> to form |u(j)>
call zdscal(h_size,beta_norm_inv,vec3,1)
!! vec3 is |u(j)>
! |w(j)> = H|u(j)>
call compute_hu(vec3,vec2,h_size)!! TODO: not implemented
!! vec2 is |w(j)>
alpha(j)=real(u_dot_v_complex(vec2,vec3,h_size))
beta(j)=beta_norm
! |v(j)> = |w(j)> - alpha(j)*|u(j)> - beta(j)*|u(j-1)>
do l=1,h_size
vec_tmp=vec2(l)-alpha(j)*vec3(l)-beta(j)*vec1(l)
vec1(l)=vec3(l)
vec3(l)=vec_tmp
enddo
!! vec1 is |u(j)>
!! vec3 is |v(j)>
enddo
end
subroutine compute_hu_hp(vec1,vec2,n_hp,h_size,spin_hp,sign_hp,idx_hp)
implicit none
integer, intent(in) :: h_size,n_hp
complex*16, intent(in) :: vec1(h_size,n_hp)
complex*16, intent(out) :: vec2(h_size,n_hp)
integer, intent(in) :: spin_hp(n_hp), idx_hp(n_hp)
double precision, intent (in) :: sign_hp(n_hp)
complex*16 :: vec1_tmp(h_size,n_hp)
integer :: i,j
BEGIN_DOC
! |vec2> = H|vec1>
!
! TODO: implement
! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}?
END_DOC
vec1_tmp(1:h_size,1:n_hp) = vec1(1:h_size,1:n_hp)
call h_u_0_hp_openmp(vec2,vec1_tmp,n_hp,h_size,spin_hp,sign_hp,idx_hp)
do j=1,n_hp
do i=1,h_size
if (cdabs(vec1_tmp(i,j) - vec1(i,j)).gt.1.d-6) then
print*,'ERROR: vec1 was changed by h_u_0_openmp'
endif
enddo
enddo
end
subroutine compute_hu(vec1,vec2,h_size)
implicit none
integer, intent(in) :: h_size
complex*16, intent(in) :: vec1(h_size)
complex*16, intent(out) :: vec2(h_size)
complex*16 :: vec1_tmp(h_size)
integer :: i
BEGIN_DOC
! |vec2> = H|vec1>
!
! TODO: implement
! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}?
END_DOC
vec1_tmp(1:h_size) = vec1(1:h_size)
call h_u_0_openmp(vec2,vec1_tmp,h_size)
do i=1,h_size
if (cdabs(vec1_tmp(i) - vec1(i)).gt.1.d-6) then
print*,'ERROR: vec1 was changed by h_u_0_openmp'
endif
enddo
end
subroutine compute_hu2(vec1,vec2,h_size)
implicit none
integer, intent(in) :: h_size
complex*16, intent(in) :: vec1(h_size)
complex*16, intent(out) :: vec2(h_size)
complex*16, allocatable :: u_tmp(:,:), s_tmp(:,:),v_tmp(:,:)
integer :: i
BEGIN_DOC
! |vec2> = H|vec1>
!
! TODO: implement
! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}?
END_DOC
allocate(u_tmp(1,h_size),s_tmp(1,h_size),v_tmp(1,h_size))
u_tmp(1,1:h_size) = vec1(1:h_size)
call h_s2_u_0_nstates_openmp(v_tmp,s_tmp,u_tmp,1,h_size)
do i=1,h_size
if (cdabs(u_tmp(1,i) - vec1(i)).gt.1.d-6) then
print*,'ERROR: vec1 was changed by h_u_0_openmp'
endif
enddo
vec2(1:h_size)=v_tmp(1,1:h_size)
deallocate(u_tmp,v_tmp,s_tmp)
end
subroutine diag_lanczos_vals_vecs(alpha, beta, nlanc, vals, vecs, sze)
implicit none
BEGIN_DOC
! diagonalization of tridiagonal form of H
! this returns eigenvalues and eigenvectors in vals,vecs
END_DOC
integer, intent(in) :: nlanc,sze
double precision, intent(in) :: alpha(sze), beta(sze)
double precision, intent(out) :: vals(sze), vecs(sze,sze)
double precision :: work(2*nlanc-2), beta_tmp(nlanc-1)
integer :: i,info
vals(1)=alpha(1)
do i=2,nlanc
vals(i)=alpha(i)
beta_tmp(i-1)=beta(i)
enddo
call dstev('V', nlanc, vals, beta_tmp, vecs, sze, work, info)
if (info.gt.0) then
print *,'WARNING: diagonalization of tridiagonal form of H did not converge'
else if (info.lt.0) then
print *,'WARNING: argument to dstev had illegal value'
endif
end
subroutine diag_lanczos_vals_hp(alpha, beta, nlanc, vals, sze,ng)
implicit none
BEGIN_DOC
! diagonalization of tridiagonal form of H
! this returns eigenvalues in vals
END_DOC
integer, intent(in) :: nlanc,sze,ng
!double precision, intent(in) :: alpha(ng,sze), beta(sze)
double precision, intent(in) :: alpha(sze,ng), beta(sze,ng)
double precision, intent(out) :: vals(sze,ng)
double precision :: work(1), beta_tmp(nlanc-1), vecs(1)
integer :: i,info,ig
do ig=1,ng
vals(1,ig)=alpha(1,ig)
do i=2,nlanc
vals(i,ig)=alpha(i,ig)
beta_tmp(i-1)=beta(i,ig)
enddo
call dstev('N', nlanc, vals(:,ig), beta_tmp, vecs, 1, work, info)
if (info.gt.0) then
print *,'WARNING: diagonalization of tridiagonal form of H did not converge'
else if (info.lt.0) then
print *,'WARNING: argument to dstev had illegal value'
endif
enddo
end
subroutine diag_lanczos_vals(alpha, beta, nlanc, vals, sze)
implicit none
BEGIN_DOC
! diagonalization of tridiagonal form of H
! this returns eigenvalues in vals
END_DOC
integer, intent(in) :: nlanc,sze
double precision, intent(in) :: alpha(sze), beta(sze)
double precision, intent(out) :: vals(sze)
double precision :: work(1), beta_tmp(nlanc-1), vecs(1)
integer :: i,info
vals(1)=alpha(1)
do i=2,nlanc
vals(i)=alpha(i)
beta_tmp(i-1)=beta(i)
enddo
call dstev('N', nlanc, vals, beta_tmp, vecs, 1, work, info)
if (info.gt.0) then
print *,'WARNING: diagonalization of tridiagonal form of H did not converge'
else if (info.lt.0) then
print *,'WARNING: argument to dstev had illegal value'
endif
end

90
src/green/plot-spec-dens.py Executable file
View File

@ -0,0 +1,90 @@
#!/bin/env python
import gzip
import sys
from math import pi
inv_pi = 1.0/pi
def spec_dens(alpha,beta,z0,g_sign,e_shift):
sze=len(alpha)
sze_b=len(beta)
if (sze != sze_b):
print('Error: size mismatch',sze,sze_b)
sys.exit(1)
z=z0-g_sign*e_shift
tmp=0.0+0.0j
#for ai,bi in zip(reversed(a),reversed(b))
for i in range(sze-1,0,-1):
tmp=-(beta[i]**2)/(z+g_sign*alpha[i]+tmp)
tmp=1.0/(z+g_sign*alpha[0]+tmp)
return -1.0 * tmp.imag * inv_pi
def printspec(ezdir,wmin,wmax,nw,eps):
gdir=ezdir+'/green/'
with open(gdir+'n_green_vec') as infile:
ngvec=int(infile.readline().strip())
with open(ezdir+'/fci/energy') as infile:
e0=float(infile.readline().strip())
with open(gdir+'n_lanczos_complete') as infile:
nlanc=int(infile.readline().strip())
with gzip.open(gdir+'green_sign.gz') as infile:
gsign0=infile.read().split()
with gzip.open(gdir+'alpha_lanczos.gz') as infile:
adata0=infile.read().split()
with gzip.open(gdir+'beta_lanczos.gz') as infile:
bdata0=infile.read().split()
adim=int(adata0.pop(0))
bdim=int(bdata0.pop(0))
gsigndim=int(gsign0.pop(0))
assert adim==2, 'dimension of alpha_lanczos should be 2'
assert bdim==2, 'dimension of beta_lanczos should be 2'
assert gsigndim==1, 'dimension of green_sign should be 1'
ngvec_2=int(gsign0.pop(0))
assert ngvec_2==ngvec, 'problem with size of green_sign.gz'
ashape=tuple(map(int,adata0[:adim]))
bshape=tuple(map(int,bdata0[:bdim]))
assert ashape==(nlanc,ngvec), 'shape of alpha_lanczos should be (nlanc, ngvec)'
assert bshape==(nlanc,ngvec), 'shape of beta_lanczos should be (nlanc, ngvec)'
amat=[]
for xi in range(ngvec):
amat.append(list(map(float,adata0[adim+xi*nlanc:adim+(xi+1)*nlanc])))
bmat=[]
b2mat=[]
for xi in range(ngvec):
#bmat.append(list(map(float,bdata0[bdim+xi*nlanc:bdim+(xi+1)*nlanc])))
b_tmp=list(map(float,bdata0[bdim+xi*nlanc:bdim+(xi+1)*nlanc]))
b2_tmp=[i*i for i in b_tmp]
bmat.append(b_tmp)
b2mat.append(b2_tmp)
gsign=list(map(float,gsign0))
dw=(wmax-wmin)/(nw-1)
wlist = [wmin+iw*dw for iw in range(nw)]
densmat=[]
for ivec in range(ngvec):
densmat.append([spec_dens(amat[ivec],bmat[ivec],iw+1.j*eps,gsign[ivec],e0) for iw in wlist])
for i,dd in enumerate(zip(*densmat)):
print(('{:15.6E}'+ngvec*'{:25.15E}').format(wlist[i],*dd))
if __name__ == '__main__':
if len(sys.argv) != 6:
print('bad args')
print('USAGE: plot-spec-dens.py ezfio omega_min omega_max n_omega epsilon')
sys.exit(1)
ezfio=sys.argv[1]
wmin=float(sys.argv[2])
wmax=float(sys.argv[3])
nw=int(sys.argv[4])
eps=float(sys.argv[5])
printspec(ezfio,wmin,wmax,nw,eps)

15
src/green/print_dets_test.irp.f Normal file
View File

@ -0,0 +1,15 @@
program print_dets_test
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
use bitmasks
implicit none
integer :: i
read*,i
print*,psi_det(:,:,i)
end

15
src/green/print_e_mo_debug.irp.f Normal file
View File

@ -0,0 +1,15 @@
program print_e_mo_debug
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
use bitmasks
implicit none
integer :: i
read*,i
call print_mo_energies(psi_det(:,:,i),N_int,mo_num)
end

178
src/green/print_h_debug.irp.f Normal file
View File

@ -0,0 +1,178 @@
program print_h_debug
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
use bitmasks
implicit none
integer :: i,j
integer, allocatable :: H_matrix_degree(:,:)
double precision, allocatable :: H_matrix_phase(:,:)
integer :: degree
integer(bit_kind), allocatable :: keys_tmp(:,:,:)
allocate(keys_tmp(N_int,2,N_det))
do i = 1, N_det
print*,''
call debug_det(psi_det(1,1,i),N_int)
do j = 1, N_int
keys_tmp(j,1,i) = psi_det(j,1,i)
keys_tmp(j,2,i) = psi_det(j,2,i)
enddo
enddo
if(N_det.gt.10000)then
print*,'Warning !!!'
print*,'Number of determinants is ',N_det
print*,'It means that the H matrix will be enormous !'
print*,'stoppping ..'
stop
endif
print*,''
print*,'Determinants '
do i = 1, N_det
enddo
allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det))
integer :: exc(0:2,2,2)
double precision :: phase
do i = 1, N_det
do j = i, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
H_matrix_degree(i,j) = degree
H_matrix_degree(j,i) = degree
phase = 0.d0
if(degree==1.or.degree==2)then
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
endif
H_matrix_phase(i,j) = phase
H_matrix_phase(j,i) = phase
enddo
enddo
print*,'H matrix '
double precision :: s2
complex*16 :: ref_h_matrix
ref_h_matrix = h_matrix_all_dets_complex(1,1)
print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion
print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion
print*,'Printing the H matrix ...'
print*,''
print*,''
!do i = 1, N_det
! H_matrix_all_dets(i,i) -= ref_h_matrix
!enddo
do i = 1, N_det
H_matrix_all_dets_complex(i,i) += nuclear_repulsion
enddo
!do i = 5,N_det
! H_matrix_all_dets(i,3) = 0.d0
! H_matrix_all_dets(3,i) = 0.d0
! H_matrix_all_dets(i,4) = 0.d0
! H_matrix_all_dets(4,i) = 0.d0
!enddo
! TODO: change for complex
do i = 1, N_det
write(*,'(I3,X,A3,2000(E24.15))')i,' | ',H_matrix_all_dets_complex(i,:)
enddo
! print*,''
! print*,''
! print*,''
! print*,'Printing the degree of excitations within the H matrix'
! print*,''
! print*,''
! do i = 1, N_det
! write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:)
! enddo
!
!
! print*,''
! print*,''
! print*,'Printing the phase of the Hamiltonian matrix elements '
! print*,''
! print*,''
! do i = 1, N_det
! write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:)
! enddo
! print*,''
! double precision, allocatable :: eigenvalues(:)
! complex*16, allocatable :: eigenvectors(:,:)
! double precision, allocatable :: s2_eigvalues(:)
! allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
! allocate (eigenvalues(N_det),s2_eigvalues(N_det))
! call lapack_diag_complex(eigenvalues,eigenvectors, &
! H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
! print*,'Two first eigenvectors '
! call u_0_S2_u_0(s2_eigvalues,eigenvectors,n_det,keys_tmp,N_int,N_det,size(eigenvectors,1))
! do j =1, N_states
! print*,'s2 = ',s2_eigvalues(j)
! print*,'e = ',eigenvalues(j)
! print*,'coefs : '
! do i = 1, N_det
! print*,'i = ',i,eigenvectors(i,j)
! enddo
! if(j>1)then
! print*,'Delta E(H) = ',eigenvalues(1) - eigenvalues(j)
! print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0
! endif
! enddo
! complex*16 :: get_mo_bielec_integral,k_a_iv,k_b_iv
! integer :: h1,p1,h2,p2
! h1 = 10
! p1 = 16
! h2 = 14
! p2 = 14
!!h1 = 1
!!p1 = 4
!!h2 = 2
!!p2 = 2
! k_a_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
! h2 = 15
! p2 = 15
! k_b_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
! print*,'k_a_iv = ',k_a_iv
! print*,'k_b_iv = ',k_b_iv
! complex*16 :: k_av,k_bv,k_ai,k_bi
! h1 = 16
! p1 = 14
! h2 = 14
! p2 = 16
! k_av = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
! h1 = 16
! p1 = 15
! h2 = 15
! p2 = 16
! k_bv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
!
! h1 = 10
! p1 = 14
! h2 = 14
! p2 = 10
! k_ai = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
!
! h1 = 10
! p1 = 15
! h2 = 15
! p2 = 10
! k_bi = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
!
! print*,'k_av, k_bv = ',k_av,k_bv
! print*,'k_ai, k_bi = ',k_ai,k_bi
! complex*16 :: k_iv
!
! h1 = 10
! p1 = 16
! h2 = 16
! p2 = 10
! k_iv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
! print*,'k_iv = ',k_iv
end

41
src/green/print_h_omp_debug.irp.f Normal file
View File

@ -0,0 +1,41 @@
program print_h_omp_debug
implicit none
read_wf = .True.
touch read_wf
call routine_omp
end
subroutine routine_omp
use bitmasks
implicit none
integer :: h_size
complex*16, allocatable :: u_tmp(:,:), s_tmp(:,:),v_tmp(:,:)
integer :: i,n_st
h_size=N_det
BEGIN_DOC
! |vec2> = H|vec1>
!
! TODO: implement
! maybe reuse parts of H_S2_u_0_nstates_{openmp,zmq}?
END_DOC
n_st=min(1000,h_size)
allocate(u_tmp(n_st,h_size),s_tmp(n_st,h_size),v_tmp(n_st,h_size))
u_tmp=(0.d0,0.d0)
v_tmp=(0.d0,0.d0)
s_tmp=(0.d0,0.d0)
do i=1,n_st
u_tmp(i,i)=(1.d0,0.d0)
enddo
call h_s2_u_0_nstates_openmp_complex(v_tmp,s_tmp,u_tmp,n_st,h_size)
do i = 1, n_st
v_tmp(i,i) += nuclear_repulsion
enddo
do i = 1, n_st
write(*,'(I3,X,A3,2000(E24.15))')i,' | ',v_tmp(i,:)
enddo
deallocate(u_tmp,v_tmp,s_tmp)
end

43
src/green/print_spectral_dens.irp.f Normal file
View File

@ -0,0 +1,43 @@
program print_spectral_dens
implicit none
BEGIN_DOC
! TODO
END_DOC
read_wf = .True.
touch read_wf
provide n_green_vec
call print_lanczos_eigvals
call print_spec
end
subroutine print_lanczos_eigvals
implicit none
integer :: i, iunit, j
integer :: getunitandopen
character(5) :: jstr
do j=1,n_green_vec
write(jstr,'(I0.3)') j
iunit = getunitandopen('lanczos_eigval_alpha_beta.out.'//trim(jstr),'w')
print *, 'printing lanczos eigenvalues, alpha, beta to "lanczos_eigval_alpha_beta.out.'//trim(jstr)//'"'
do i=1,n_lanczos_iter
write(iunit,'(I6,3(E25.15))') i, lanczos_eigvals(i,j), alpha_lanczos(i,j), beta_lanczos(i,j)
enddo
close(iunit)
enddo
end
subroutine print_spec
implicit none
integer :: i, iunit, j
integer :: getunitandopen
character(5) :: jstr
do j=1,n_green_vec
write(jstr,'(I0.3)') j
iunit = getunitandopen('omega_A.out.'//trim(jstr),'w')
print *, 'printing frequency, spectral density to "omega_A.out.'//trim(jstr)//'"'
do i=1,n_omega
write(iunit,'(2(E25.15))') omega_list(i), spectral_lanczos(i,j)
enddo
close(iunit)
enddo
end

614
src/green/utils_hp.irp.f Normal file
View File

@ -0,0 +1,614 @@
subroutine print_mo_energies(key_ref,nint,nmo)
use bitmasks
BEGIN_DOC
! get mo energies for one det
END_DOC
implicit none
integer, intent(in) :: nint, nmo
integer(bit_kind), intent(in) :: key_ref(nint,2)
double precision, allocatable :: e_mo(:,:)
integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2)
integer :: n_occ(2), n_virt(2)
integer, parameter :: int_spin2(1:2) = (/2,1/)
integer :: i,j,ispin,jspin,i0,j0,k
integer(bit_kind), allocatable :: key_virt(:,:)
integer, allocatable :: is_occ(:,:)
allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2),key_virt(nint,2),e_mo(nmo,2),is_occ(nmo,2))
is_occ=0
call bitstring_to_list_ab(key_ref,occ,n_occ,nint)
do i=1,nint
do ispin=1,2
key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin))
enddo
enddo
call bitstring_to_list_ab(key_virt,virt,n_virt,nint)
e_mo(1:nmo,1)=mo_one_e_integrals_diag(1:nmo)
e_mo(1:nmo,2)=mo_one_e_integrals_diag(1:nmo)
do ispin=1,2
jspin=int_spin2(ispin)
do i0=1,n_occ(ispin)
i=occ(i0,ispin)
is_occ(i,ispin)=1
do j0=i0+1,n_occ(ispin)
j=occ(j0,ispin)
e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj_anti(i,j)
e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j)
enddo
do k=2,ispin
do j0=1,n_occ(jspin)
j=occ(j0,jspin)
e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj(i,j)
e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) !can delete this and remove k level of loop
enddo
enddo
do j0=1,n_virt(ispin)
j=virt(j0,ispin)
e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j)
enddo
do j0=1,n_virt(jspin)
j=virt(j0,jspin)
e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j)
enddo
enddo
enddo
do i=1,nmo
write(6,'(2(I5),2(E25.15))')is_occ(i,1),is_occ(i,2),e_mo(i,1),e_mo(i,2)
enddo
deallocate(occ,virt,key_virt,e_mo,is_occ)
end
subroutine get_mo_energies(key_ref,nint,nmo,e_mo)
use bitmasks
BEGIN_DOC
! get mo energies for one det
END_DOC
implicit none
integer, intent(in) :: nint, nmo
integer(bit_kind), intent(in) :: key_ref(nint,2)
double precision, intent(out) :: e_mo(nmo,2)
integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2)
integer :: n_occ(2), n_virt(2)
integer, parameter :: int_spin2(1:2) = (/2,1/)
integer :: i,j,ispin,jspin,i0,j0,k
integer(bit_kind), allocatable :: key_virt(:,:)
allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2),key_virt(nint,2))
call bitstring_to_list_ab(key_ref,occ,n_occ,nint)
do i=1,nint
do ispin=1,2
key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin))
enddo
enddo
call bitstring_to_list_ab(key_virt,virt,n_virt,nint)
e_mo(1:nmo,1)=mo_one_e_integrals_diag(1:nmo)
e_mo(1:nmo,2)=mo_one_e_integrals_diag(1:nmo)
do ispin=1,2
jspin=int_spin2(ispin)
do i0=1,n_occ(ispin)
i=occ(i0,ispin)
do j0=i0+1,n_occ(ispin)
j=occ(j0,ispin)
e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj_anti(i,j)
e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j)
enddo
do k=2,ispin
do j0=1,n_occ(jspin)
j=occ(j0,jspin)
e_mo(i,ispin) = e_mo(i,ispin) + mo_two_e_integrals_jj(i,j)
e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j) !can delete this and remove k level of loop
enddo
enddo
do j0=1,n_virt(ispin)
j=virt(j0,ispin)
e_mo(j,ispin) = e_mo(j,ispin) + mo_two_e_integrals_jj_anti(i,j)
enddo
do j0=1,n_virt(jspin)
j=virt(j0,jspin)
e_mo(j,jspin) = e_mo(j,jspin) + mo_two_e_integrals_jj(i,j)
enddo
enddo
enddo
deallocate(occ,virt,key_virt)
end
subroutine get_mask_phase_new(det1, pm, Nint)
use bitmasks
BEGIN_DOC
! phasemask copied from qp2
! return phasemask of det1 in pm
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2)
integer(bit_kind), intent(out) :: pm(Nint,2)
integer(bit_kind) :: tmp1, tmp2
integer :: i
pm(1:Nint,1:2) = det1(1:Nint,1:2)
tmp1 = 0_8
tmp2 = 0_8
do i=1,Nint
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1))
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2))
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4))
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8))
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16))
pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32))
pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32))
pm(i,1) = ieor(pm(i,1), tmp1)
pm(i,2) = ieor(pm(i,2), tmp2)
if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1)
if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2)
end do
end subroutine
subroutine get_phase_hp(g_idx_int,g_idx_bit,g_spin,g_sign,det_in,g_det_phase,nint,n_g)
use bitmasks
implicit none
integer, intent(in) :: nint,n_g
integer, intent(in) :: g_idx_int(n_g), g_idx_bit(n_g),g_spin(n_g)
double precision, intent(in) :: g_sign(n_g)
integer(bit_kind), intent(in) :: det_in(nint,2)
double precision, intent(out) :: g_det_phase(n_g)
integer(bit_kind) :: tmp_spindet(nint), pm(nint,2)
double precision, parameter :: phase_dble(0:1) = (/1.d0,-1.d0/)
integer :: i
logical :: is_allowed(n_g), all_banned, is_filled
all_banned=.True.
do i=1,n_g
tmp_spindet(1:nint) = det_in(1:nint,g_spin(i))
call spinorb_is_filled_int_bit(tmp_spindet,g_idx_int(i),g_idx_bit(i),nint,is_filled)
is_allowed(i) = (.not.(((g_sign(i)<0).and.(.not.is_filled)).or.((g_sign(i)>0).and.(is_filled))))
all_banned=(all_banned.and.(.not.is_allowed(i)))
enddo
if (all_banned) then
g_det_phase(:)=0.d0
else
call get_mask_phase_new(det_in,pm,nint)
do i=1,n_g
if (is_allowed(i)) then
g_det_phase(i) = phase_dble(popcnt(iand(ibset(0_bit_kind,g_idx_bit(i)),pm(g_idx_int(i),g_spin(i)))))
else
g_det_phase(i)=0.d0
endif
enddo
endif
end
subroutine get_homo_lumo(key_ref,nint,nmo,idx_homo_lumo,spin_homo_lumo)
use bitmasks
implicit none
integer, intent(in) :: nint,nmo
integer(bit_kind), intent(in) :: key_ref(nint,2)
integer, intent(out) :: idx_homo_lumo(2), spin_homo_lumo(2)
double precision, allocatable :: e_mo(:,:)
integer, allocatable :: occ(:,:),virt(:,:) !(nint*bit_kind_size,2)
integer :: n_occ(2), n_virt(2)
integer :: i,i0,ispin
integer(bit_kind), allocatable :: key_virt(:,:)
double precision :: maxocc(2), minvirt(2)
integer :: imaxocc(2), iminvirt(2)
allocate(e_mo(nmo,2),key_virt(nint,2),occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2))
call get_mo_energies(key_ref,nint,nmo,e_mo)
!allocate(occ(nint*bit_kind_size,2),virt(nint*bit_kind_size,2))
call bitstring_to_list_ab(key_ref,occ,n_occ,nint)
do i=1,nint
do ispin=1,2
key_virt(i,ispin)=xor(full_ijkl_bitmask(i),key_ref(i,ispin))
enddo
enddo
call bitstring_to_list_ab(key_virt,virt,n_virt,nint)
maxocc=-1.d20 !maybe use -1.d0*huge(0.d0)?
minvirt=1.d20
imaxocc=-1
iminvirt=-1
do ispin=1,2
do i0=1,n_occ(ispin)
i=occ(i0,ispin)
if (e_mo(i,ispin).gt.maxocc(ispin)) then
maxocc(ispin)=e_mo(i,ispin)
imaxocc(ispin)=i
endif
enddo
do i0=1,n_virt(ispin)
i=virt(i0,ispin)
if (e_mo(i,ispin).lt.minvirt(ispin)) then
minvirt(ispin)=e_mo(i,ispin)
iminvirt(ispin)=i
endif
enddo
enddo
double precision :: e_mo_thresh
e_mo_thresh = 1.d-8
!these should both just be 2x2 arrays, but performance here doesn't really matter and this is more readable
!if (maxocc(1).ge.maxocc(2)) then
if ((maxocc(2)-maxocc(1)).le.e_mo_thresh) then
spin_homo_lumo(1)=1
else
spin_homo_lumo(1)=2
endif
if ((minvirt(1)-minvirt(2)).le.e_mo_thresh) then
spin_homo_lumo(2)=1
else
spin_homo_lumo(2)=2
endif
idx_homo_lumo(1)=imaxocc(spin_homo_lumo(1))
idx_homo_lumo(2)=iminvirt(spin_homo_lumo(2))
deallocate(e_mo,occ,virt,key_virt)
end
subroutine get_list_hp_banned_ab(tmp_det,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,nint,all_banned)
use bitmasks
implicit none
BEGIN_DOC
! input determinant tmp_det and list of single holes/particles
! for each hole/particle, determine whether it is filled/empty in tmp_det
! return which are disallowed in exc_is_banned
! if all are banned, set all_banned to true
END_DOC
integer, intent(in) :: N_hp,nint
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
integer(bit_kind), intent(in) :: tmp_det(nint,2)
logical, intent(out) :: exc_is_banned(N_hp)
logical, intent(out) :: all_banned
integer :: i
logical :: is_filled
all_banned = .True.
do i=1,N_hp
call orb_is_filled(tmp_det,idx_hp(i),spin_hp(i),nint,is_filled)
if (sign_hp(i).gt.0) then ! particle creation, banned if already filled
exc_is_banned(i) = is_filled
else ! hole creation, banned if already empty
exc_is_banned(i) = (.not.is_filled)
endif
all_banned = (all_banned.and.exc_is_banned(i))
enddo
end
subroutine get_list_hp_banned_single_spin(tmp_spindet,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,ispin,nint,all_banned)
use bitmasks
implicit none
BEGIN_DOC
! input spindeterminant tmp_spindet and list of single holes/particles
! tmp_spindet is only one spin part of a full det, with spin==ispin
! for each hole/particle, determine whether it is filled/empty in tmp_det
! return which are disallowed in exc_is_banned
! if all are banned, set all_banned to true
END_DOC
integer, intent(in) :: N_hp, ispin, nint
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
integer(bit_kind), intent(in) :: tmp_spindet(nint)
logical, intent(out) :: exc_is_banned(N_hp)
logical, intent(out) :: all_banned
integer :: i
logical :: is_filled
all_banned = .True.
do i=1,N_hp
if (spin_hp(i).eq.ispin) then
call orb_is_filled_single_spin(tmp_spindet,idx_hp(i),nint,is_filled)
if (sign_hp(i).gt.0) then ! particle creation, banned if already filled
exc_is_banned(i) = is_filled
else ! hole creation, banned if already empty
exc_is_banned(i) = (.not.is_filled)
endif
else
exc_is_banned(i) = .False.
endif
all_banned = (all_banned.and.exc_is_banned(i))
enddo
end
subroutine get_list_hp_banned_spin(tmp_det,N_hp,exc_is_banned,spin_hp,sign_hp,idx_hp,ispin,nint,all_banned)
use bitmasks
implicit none
BEGIN_DOC
! input determinant tmp_det and list of single holes/particles
! for each hole/particle, determine whether it is filled/empty in tmp_det
! return which are disallowed in exc_is_banned
! if all are banned, set all_banned to true
! only consider tmp_det(1:N_int, ispin)
END_DOC
integer, intent(in) :: N_hp, ispin, nint
integer, intent(in) :: spin_hp(N_hp), idx_hp(N_hp)
double precision, intent(in) :: sign_hp(N_hp)
integer(bit_kind), intent(in) :: tmp_det(nint,2)
logical, intent(out) :: exc_is_banned(N_hp)
logical, intent(out) :: all_banned
integer(bit_kind) :: spindet(nint)
integer :: i
logical :: is_filled
spindet(1:nint) = tmp_det(1:nint,ispin)
all_banned = .True.
do i=1,N_hp
if (spin_hp(i).eq.ispin) then
call orb_is_filled(tmp_det,idx_hp(i),ispin,nint,is_filled)
if (sign_hp(i).gt.0) then ! particle creation, banned if already filled
exc_is_banned(i) = is_filled
else ! hole creation, banned if already empty
exc_is_banned(i) = (.not.is_filled)
endif
else
exc_is_banned(i) = .False.
endif
all_banned = (all_banned.and.exc_is_banned(i))
enddo
end
subroutine spinorb_is_filled_int_bit(key_ref,iorb_int,iorb_bit,Nint,is_filled)
use bitmasks
implicit none
BEGIN_DOC
! determine whether iorb is filled in key_ref
! iorb is specified by int and bit locations within the determinant
END_DOC
integer, intent(in) :: iorb_int, iorb_bit, Nint
integer(bit_kind), intent(in) :: key_ref(Nint)
logical, intent(out) :: is_filled
ASSERT (iorb_int > 0)
ASSERT (iorb_bit >= 0)
ASSERT (Nint > 0)
is_filled = btest(key_ref(iorb_int),iorb_bit)
end
subroutine orb_is_filled_int_bit(key_ref,iorb_int,iorb_bit,ispin,Nint,is_filled)
use bitmasks
implicit none
BEGIN_DOC
! todo: not finished
! determine whether iorb is filled in key_ref
! iorb is specified by int and bit locations within the determinant
END_DOC
integer, intent(in) :: iorb_int, iorb_bit, ispin, Nint
integer(bit_kind), intent(in) :: key_ref(Nint,2)
logical, intent(out) :: is_filled
ASSERT (iorb_int > 0)
ASSERT (iorb_bit >= 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
is_filled = btest(key_ref(iorb_int,ispin),iorb_bit)
! call spinorb_is_filled_int_bit(key_ref(1,ispin),iorb_int,iorb_bit,Nint,is_filled)
end
subroutine get_orb_int_bit(iorb,iorb_int,iorb_bit)
BEGIN_DOC
! get int and bit corresponding to orbital index iorb
END_DOC
use bitmasks
implicit none
integer, intent(in) :: iorb
integer, intent(out) :: iorb_int, iorb_bit
ASSERT (iorb > 0)
iorb_int = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (iorb_int > 0)
iorb_bit = iorb - ishft(iorb_int-1,bit_kind_shift)-1
ASSERT (iorb_bit >= 0)
end
subroutine orb_is_filled_single_spin(key_ref,iorb,Nint,is_filled)
use bitmasks
implicit none
BEGIN_DOC
! determine whether iorb is filled in key_ref
! key_ref is single alpha or beta determinant
END_DOC
integer, intent(in) :: iorb, Nint
integer(bit_kind), intent(in) :: key_ref(Nint)
logical, intent(out) :: is_filled
integer :: k,l
ASSERT (iorb > 0)
ASSERT (Nint > 0)
! k is index of the int where iorb is found
! l is index of the bit where iorb is found
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k >0)
l = iorb - ishft(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
is_filled = btest(key_ref(k),l)
end
subroutine orb_is_filled(key_ref,iorb,ispin,Nint,is_filled)
use bitmasks
implicit none
BEGIN_DOC
! determine whether iorb, ispin is filled in key_ref
! key_ref has alpha and beta parts
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer(bit_kind), intent(in) :: key_ref(Nint,2)
logical, intent(out) :: is_filled
integer :: k,l
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
! k is index of the int where iorb is found
! l is index of the bit where iorb is found
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k >0)
l = iorb - ishft(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
is_filled = btest(key_ref(k,ispin),l)
end
subroutine ac_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase)
use bitmasks
implicit none
BEGIN_DOC
! apply creation operator to key_ref
! add electron with spin ispin to orbital with index iorb
! output resulting det and phase in key_new and phase
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer(bit_kind), intent(in) :: key_ref(Nint,2)
integer(bit_kind), intent(out) :: key_new(Nint,2)
double precision, intent(out) :: phase
integer :: k,l,i
double precision, parameter :: p(0:1) = (/ 1.d0, -1.d0 /)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
key_new=key_ref
! alpha det is list of Nint 64-bit ints
! k is index of the int where iorb is found
! l is index of the bit where iorb is found
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k >0)
l = iorb - ishft(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
key_new(k,ispin) = ibset(key_new(k,ispin),l)
integer(bit_kind) :: parity_filled
! I assume here that the ordering is all alpha spinorbs and then all beta spinorbs
! If we add an alpha electron, parity is not affected by beta part of determinant
! (only need number of alpha occupied orbs below iorb)
! If we add a beta electron, the parity is affected by alpha part
! (need total number of occupied alpha orbs (all of which come before beta)
! and total number of beta occupied orbs below iorb)
if (ispin==1) then
parity_filled=0_bit_kind
else
parity_filled=iand(int(elec_alpha_num,bit_kind),1_bit_kind)
endif
! get parity due to orbs in other ints (with lower indices)
do i=1,k-1
parity_filled = iand(int(popcnt(key_ref(i,ispin)),bit_kind),parity_filled)
enddo
! get parity due to orbs in same int as iorb
! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0
parity_filled = iand(int(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),bit_kind),parity_filled)
phase = p(iand(1_bit_kind,parity_filled))
end
subroutine a_operator_phase(key_new,key_ref,iorb,ispin,Nint,phase)
use bitmasks
implicit none
BEGIN_DOC
! apply annihilation operator to key_ref
! remove electron with spin ispin to orbital with index iorb
! output resulting det and phase in key_new and phase
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer(bit_kind), intent(in) :: key_ref(Nint,2)
integer(bit_kind), intent(out) :: key_new(Nint,2)
double precision, intent(out) :: phase
integer :: k,l,i
double precision, parameter :: p(0:1) = (/ 1.d0, -1.d0 /)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
key_new=key_ref
! alpha det is list of Nint 64-bit ints
! k is index of the int where iorb is found
! l is index of the bit where iorb is found
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k >0)
l = iorb - ishft(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
key_new(k,ispin) = ibclr(key_new(k,ispin),l)
integer(bit_kind) :: parity_filled
! I assume here that the ordering is all alpha spinorbs and then all beta spinorbs
! If we add an alpha electron, parity is not affected by beta part of determinant
! (only need number of alpha occupied orbs below iorb)
! If we add a beta electron, the parity is affected by alpha part
! (need total number of occupied alpha orbs (all of which come before beta)
! and total number of beta occupied orbs below iorb)
if (ispin==1) then
parity_filled=0_bit_kind
else
parity_filled=iand(int(elec_alpha_num,bit_kind),1_bit_kind)
endif
! get parity due to orbs in other ints (with lower indices)
do i=1,k-1
parity_filled = iand(int(popcnt(key_ref(i,ispin)),bit_kind),parity_filled)
enddo
! get parity due to orbs in same int as iorb
! ishft(1_bit_kind,l)-1 has its l rightmost bits set to 1, other bits set to 0
parity_filled = iand(int(popcnt(iand(ishft(1_bit_kind,l)-1,key_ref(k,ispin))),bit_kind),parity_filled)
phase = p(iand(1_bit_kind,parity_filled))
end
!BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_diag,(mo_num)]
! implicit none
! integer :: i
! BEGIN_DOC
! ! diagonal elements of mo_mono_elec_integral array
! END_DOC
! print*,'Providing the mono electronic integrals (diagonal)'
!
! do i = 1, mo_num
! mo_mono_elec_integral_diag(i) = real(mo_mono_elec_integral(i,i))
! enddo
!
!END_PROVIDER

6
src/mo_two_e_ints/df_mo_ints.irp.f
View File

@ -48,7 +48,8 @@ subroutine mo_map_fill_from_df_dot
logical :: use_map1
integer(key_kind) :: idx_tmp
double precision :: sign
complex*16, external :: zdotc
!complex*16, external :: zdotc
complex*16, external :: zdotu
mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt
@ -145,7 +146,8 @@ subroutine mo_map_fill_from_df_dot
if ((j==l) .and. (i>k)) exit
call idx2_tri_int(i,k,ik2)
if (ik2 > jl2) exit
integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
!integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
integral = zdotu(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
! print*,i,k,j,l,real(integral),imag(integral)
if (cdabs(integral) < mo_integrals_threshold) then
cycle

91
src/tools/fcidump.irp.f
View File

@ -18,6 +18,97 @@ program fcidump
! electrons
!
END_DOC
if (is_complex) then
call fcidump_complex
else
call fcidump_real
endif
end
subroutine fcidump_complex
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.FCIDUMP'
i_unit_output = getUnitAndOpen(output,'w')
integer :: i,j,k,l
integer :: i1,j1,k1,l1
integer :: i2,j2,k2,l2,ik2,jl2
integer :: ki,kj,kk,kl
integer :: ii,ij,ik,il
integer*8 :: m
character*(2), allocatable :: A(:)
write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, &
', MS2=', (elec_alpha_num-elec_beta_num), ','
allocate (A(n_act_orb))
A = '1,'
write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb)
write(i_unit_output,*) 'ISYM=0,'
write(i_unit_output,*) '/'
deallocate(A)
integer(key_kind), allocatable :: keys(:)
double precision, allocatable :: values(:)
integer(cache_map_size_kind) :: n_elements, n_elements_max
PROVIDE mo_two_e_integrals_in_map
complex*16 :: get_two_e_integral_complex, integral
do kl=1,kpt_num
do kj=1,kl
do kk=1,kl
ki=kconserv(kl,kk,kj)
if (ki>kl) cycle
do l1=1,n_act_orb_kpts(kl)
il=list_act_kpts(l1,kl)
l = (kl-1)*mo_num_per_kpt + il
do j1=1,n_act_orb_kpts(kj)
ij=list_act_kpts(j1,kj)
j = (kj-1)*mo_num_per_kpt + ij
if (j>l) exit
call idx2_tri_int(j,l,jl2)
do k1=1,n_act_orb_kpts(kk)
ik=list_act_kpts(k1,kk)
k = (kk-1)*mo_num_per_kpt + ik
if (k>l) exit
do i1=1,n_act_orb_kpts(ki)
ii=list_act_kpts(i1,ki)
i = (ki-1)*mo_num_per_kpt + ii
if ((j==l) .and. (i>k)) exit
call idx2_tri_int(i,k,ik2)
if (ik2 > jl2) exit
integral = get_two_e_integral_complex(i,j,k,l,mo_integrals_map,mo_integrals_map_2)
if (cdabs(integral) > mo_integrals_threshold) then
write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral), dimag(integral),i,k,j,l
endif
enddo
enddo
enddo
enddo
enddo
enddo
enddo
do kj=1,kpt_num
do j1=1,n_act_orb_kpts(kj)
ij = list_act_kpts(j1,kj)
j = (kj-1)*mo_num_per_kpt + ij
do i1=j1,n_act_orb_kpts(kj)
ii = list_act_kpts(i1,kj)
i = (kj-1)*mo_num_per_kpt + ii
integral = mo_one_e_integrals_kpts(ii,ij,kj) + core_fock_operator_complex(i,j)
if (cdabs(integral) > mo_integrals_threshold) then
write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral),dimag(integral), i,j,0,0
endif
enddo
enddo
enddo
write(i_unit_output,*) core_energy, 0, 0, 0, 0
end
subroutine fcidump_real
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.FCIDUMP'

1
src/utils/constants.include.F
View File

@ -7,6 +7,7 @@ double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))
double precision, parameter :: pi_5_2 = 34.9868366552d0
double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0)
double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
double precision, parameter :: inv_pi = 1.d0/dacos(-1.d0)
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
double precision, parameter :: thresh = 1.d-15