10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-22 13:12:23 +02:00
QuantumPackage/plugins/local/bi_ort_ints/bi_ort_ints.irp.f

476 lines
11 KiB
Fortran
Raw Permalink Normal View History

2023-06-04 09:58:29 +02:00
! ---
2023-02-07 13:27:19 +01:00
program bi_ort_ints
2023-06-04 09:58:29 +02:00
2023-02-07 13:27:19 +01:00
BEGIN_DOC
2023-06-04 09:58:29 +02:00
! TODO : Put the documentation of the program here
2023-02-07 13:27:19 +01:00
END_DOC
2023-06-04 09:58:29 +02:00
implicit none
2023-02-07 13:27:19 +01:00
my_grid_becke = .True.
2023-07-02 21:49:25 +02:00
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
2023-06-04 09:58:29 +02:00
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
2023-06-02 09:11:32 +02:00
! call test_3e
2023-06-03 22:12:30 +02:00
! call test_5idx
! call test_5idx2
call test_4idx()
2023-09-06 21:03:22 +02:00
!call test_4idx_n4()
!call test_4idx2()
!call test_5idx2
!call test_5idx
2023-09-06 21:03:22 +02:00
2023-06-02 11:40:21 +02:00
end
subroutine test_5idx2
PROVIDE three_e_5_idx_cycle_2_bi_ort
2023-02-07 13:27:19 +01:00
end
2023-06-04 09:58:29 +02:00
subroutine test_4idx2()
2023-06-04 15:27:07 +02:00
!PROVIDE three_e_4_idx_direct_bi_ort
PROVIDE three_e_4_idx_exch23_bi_ort
2023-06-04 09:58:29 +02:00
end
2023-02-07 13:27:19 +01:00
subroutine test_3e
implicit none
integer :: i,k,j,l,m,n,ipoint
double precision :: accu, contrib,new,ref
i = 1
k = 1
2023-06-02 00:33:37 +02:00
n = 0
2023-02-07 13:27:19 +01:00
accu = 0.d0
do i = 1, mo_num
2023-06-02 09:11:32 +02:00
do k = 1, mo_num
2023-02-07 13:27:19 +01:00
do j = 1, mo_num
2023-06-02 09:11:32 +02:00
do l = 1, mo_num
2023-02-07 13:27:19 +01:00
do m = 1, mo_num
2023-06-02 09:11:32 +02:00
do n = 1, mo_num
call give_integrals_3_body_bi_ort(n, l, k, m, j, i, new)
call give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, ref)
2023-02-07 13:27:19 +01:00
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'pb !!'
print*,i,k,j,l,m,n
print*,ref,new,contrib
2023-06-02 00:33:37 +02:00
stop
2023-02-07 13:27:19 +01:00
endif
2023-06-02 09:11:32 +02:00
enddo
2023-02-07 13:27:19 +01:00
enddo
enddo
enddo
enddo
enddo
print*,'accu = ',accu/dble(mo_num)**6
end
2023-06-02 09:11:32 +02:00
subroutine test_5idx
implicit none
integer :: i,k,j,l,m,n,ipoint
double precision :: accu, contrib,new,ref
2023-06-02 20:32:31 +02:00
double precision, external :: three_e_5_idx_exch12_bi_ort
2023-06-02 09:11:32 +02:00
i = 1
k = 1
n = 0
accu = 0.d0
2023-06-04 16:45:38 +02:00
PROVIDE three_e_5_idx_direct_bi_ort_old
2023-06-02 09:11:32 +02:00
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
2023-06-02 20:32:31 +02:00
! if (dabs(three_e_5_idx_direct_bi_ort(m,l,j,k,i) - three_e_5_idx_exch12_bi_ort(m,l,i,k,j)) > 1.d-10) then
! stop
! endif
2023-06-04 16:45:38 +02:00
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'direct'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
!
! new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
2023-06-02 20:32:31 +02:00
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
2023-06-04 16:45:38 +02:00
! print*,'exch12'
2023-06-02 20:32:31 +02:00
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
!
! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
2023-06-02 10:34:05 +02:00
2023-06-02 20:32:31 +02:00
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle2'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch23'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch13'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle2'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch23'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch13'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
2023-06-02 09:11:32 +02:00
enddo
enddo
enddo
enddo
enddo
print*,'accu = ',accu/dble(mo_num)**5
end
2023-06-03 22:12:30 +02:00
! ---
subroutine test_4idx_n4()
implicit none
integer :: i, j, k, l
double precision :: accu, contrib, new, ref, thr
thr = 1d-10
PROVIDE three_e_4_idx_direct_bi_ort_old
PROVIDE three_e_4_idx_direct_bi_ort_n4
accu = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_direct_bi_ort_n4 (l,k,j,i)
ref = three_e_4_idx_direct_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_direct_bi_ort_n4'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_direct_bi_ort_n4 = ', accu / dble(mo_num)**4
! ---
PROVIDE three_e_4_idx_exch13_bi_ort_old
PROVIDE three_e_4_idx_exch13_bi_ort_n4
accu = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_exch13_bi_ort_n4 (l,k,j,i)
ref = three_e_4_idx_exch13_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch13_bi_ort_n4'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_exch13_bi_ort_n4 = ', accu / dble(mo_num)**4
! ---
PROVIDE three_e_4_idx_cycle_1_bi_ort_old
PROVIDE three_e_4_idx_cycle_1_bi_ort_n4
accu = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_cycle_1_bi_ort_n4 (l,k,j,i)
ref = three_e_4_idx_cycle_1_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_cycle_1_bi_ort_n4'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_cycle_1_bi_ort_n4 = ', accu / dble(mo_num)**4
! ---
PROVIDE three_e_4_idx_exch23_bi_ort_old
PROVIDE three_e_4_idx_exch23_bi_ort_n4
accu = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_exch23_bi_ort_n4 (l,k,j,i)
ref = three_e_4_idx_exch23_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch23_bi_ort_n4'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_exch23_bi_ort_n4 = ', accu / dble(mo_num)**4
! ---
return
end
! ---
2023-06-03 22:12:30 +02:00
subroutine test_4idx()
implicit none
integer :: i, j, k, l
2023-09-06 21:03:22 +02:00
double precision :: accu, contrib, new, ref, thr, norm
2023-06-03 22:12:30 +02:00
thr = 1d-10
2023-06-03 22:12:30 +02:00
PROVIDE three_e_4_idx_direct_bi_ort_old
PROVIDE three_e_4_idx_direct_bi_ort
accu = 0.d0
2023-09-06 21:03:22 +02:00
norm = 0.d0
2023-06-03 22:12:30 +02:00
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_direct_bi_ort (l,k,j,i)
ref = three_e_4_idx_direct_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_direct_bi_ort'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
2023-09-06 21:03:22 +02:00
accu += contrib
norm += dabs(ref)
2023-06-03 22:12:30 +02:00
enddo
enddo
enddo
enddo
2023-09-06 21:03:22 +02:00
print*, ' accu on three_e_4_idx_direct_bi_ort (%) = ', 100.d0 * accu / norm
2023-06-03 22:12:30 +02:00
! ---
PROVIDE three_e_4_idx_exch13_bi_ort_old
PROVIDE three_e_4_idx_exch13_bi_ort
accu = 0.d0
2023-09-06 21:03:22 +02:00
norm = 0.d0
2023-06-03 22:12:30 +02:00
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_exch13_bi_ort (l,k,j,i)
ref = three_e_4_idx_exch13_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch13_bi_ort'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
2023-09-06 21:03:22 +02:00
accu += contrib
norm += dabs(ref)
2023-06-03 22:12:30 +02:00
enddo
enddo
enddo
enddo
2023-09-06 21:03:22 +02:00
print*, ' accu on three_e_4_idx_exch13_bi_ort (%) = ', 100.d0 * accu / norm
2023-06-03 22:12:30 +02:00
! ---
PROVIDE three_e_4_idx_cycle_1_bi_ort_old
PROVIDE three_e_4_idx_cycle_1_bi_ort
accu = 0.d0
2023-09-06 21:03:22 +02:00
norm = 0.d0
2023-06-03 22:12:30 +02:00
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_cycle_1_bi_ort (l,k,j,i)
ref = three_e_4_idx_cycle_1_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_cycle_1_bi_ort'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
2023-09-06 21:03:22 +02:00
accu += contrib
norm += dabs(ref)
2023-06-03 22:12:30 +02:00
enddo
enddo
enddo
enddo
2023-09-06 21:03:22 +02:00
print*, ' accu on three_e_4_idx_cycle_1_bi_ort (%) = ', 100.d0 * accu / norm
2023-06-03 22:12:30 +02:00
! ---
PROVIDE three_e_4_idx_exch23_bi_ort_old
PROVIDE three_e_4_idx_exch23_bi_ort
accu = 0.d0
2023-09-06 21:03:22 +02:00
norm = 0.d0
2023-06-03 22:12:30 +02:00
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
new = three_e_4_idx_exch23_bi_ort (l,k,j,i)
ref = three_e_4_idx_exch23_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch23_bi_ort'
print*, l, k, j, i
print*, ref, new, contrib
stop
endif
2023-09-06 21:03:22 +02:00
accu += contrib
norm += dabs(ref)
2023-06-03 22:12:30 +02:00
enddo
enddo
enddo
enddo
2023-09-06 21:03:22 +02:00
print*, ' accu on three_e_4_idx_exch23_bi_ort (%) = ', 100.d0 * accu / norm
2023-06-03 22:12:30 +02:00
! ---
return
end
2023-09-06 21:03:22 +02:00