mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-23 03:53:29 +01:00
rm deb code from HF module
This commit is contained in:
AbdAmmar
parent
398ca5ceb7
commit
9071a64ce1
@ -1,909 +0,0 @@
|
|||||||
|
|
||||||
program deb_ao_2e_int
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
!call main()
|
|
||||||
call check_ao_one_e_integral_cgtos()
|
|
||||||
!call check_ao_two_e_integral_cgtos()
|
|
||||||
!call check_crint1()
|
|
||||||
!call check_crint2()
|
|
||||||
!call check_crint3()
|
|
||||||
!call check_crint4()
|
|
||||||
!call check_crint5()
|
|
||||||
!call check_crint6()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine main()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i, j
|
|
||||||
|
|
||||||
|
|
||||||
PROVIDE ao_overlap
|
|
||||||
PROVIDE ao_kinetic_integrals
|
|
||||||
PROVIDE ao_integrals_n_e
|
|
||||||
|
|
||||||
print*, "ao_overlap:"
|
|
||||||
do i = 1, ao_num
|
|
||||||
print*, (ao_overlap(i,j), j=1, ao_num)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "ao_kinetic_integrals:"
|
|
||||||
do i = 1, ao_num
|
|
||||||
print*, (ao_kinetic_integrals(i,j), j=1, ao_num)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "ao_integrals_n_e:"
|
|
||||||
do i = 1, ao_num
|
|
||||||
print*, (ao_integrals_n_e(i,j), j=1, ao_num)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
return
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_ao_one_e_integral_cgtos()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i, j
|
|
||||||
double precision :: acc, nrm, dif
|
|
||||||
double precision :: tmp1, tmp2
|
|
||||||
double precision :: t1, t2, tt
|
|
||||||
|
|
||||||
PROVIDE ao_overlap ao_overlap_cgtos
|
|
||||||
PROVIDE ao_integrals_n_e ao_integrals_n_e_cgtos
|
|
||||||
PROVIDE ao_kinetic_integrals ao_kinetic_integrals_cgtos
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
! print *, "overlap:"
|
|
||||||
! acc = 0.d0
|
|
||||||
! nrm = 0.d0
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! tmp1 = ao_overlap (i,j)
|
|
||||||
! tmp2 = ao_overlap_cgtos(i,j)
|
|
||||||
! dif = abs(tmp1 - tmp2)
|
|
||||||
! if(dif .gt. 1d-10) then
|
|
||||||
! print*, ' error on:', i, j
|
|
||||||
! print*, tmp1, tmp2, dif
|
|
||||||
! !stop
|
|
||||||
! endif
|
|
||||||
! acc += dif
|
|
||||||
! nrm += abs(tmp1)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! print *, ' acc (%) = ', 100.d0 * acc / nrm
|
|
||||||
!
|
|
||||||
! ! ---
|
|
||||||
!
|
|
||||||
! print *, "kinetic:"
|
|
||||||
! acc = 0.d0
|
|
||||||
! nrm = 0.d0
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! tmp1 = ao_kinetic_integrals (i,j)
|
|
||||||
! tmp2 = ao_kinetic_integrals_cgtos(i,j)
|
|
||||||
! dif = abs(tmp1 - tmp2)
|
|
||||||
! if(dif .gt. 1d-10) then
|
|
||||||
! print*, ' error on:', i, j
|
|
||||||
! print*, tmp1, tmp2, dif
|
|
||||||
! !stop
|
|
||||||
! endif
|
|
||||||
! acc += dif
|
|
||||||
! nrm += abs(tmp1)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! print *, ' acc (%) = ', 100.d0 * acc / nrm
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
print *, "NAI:"
|
|
||||||
acc = 0.d0
|
|
||||||
nrm = 0.d0
|
|
||||||
do i = 1, ao_num
|
|
||||||
!do i = 9, 9
|
|
||||||
do j = 1, ao_num
|
|
||||||
!do j = 16, 16
|
|
||||||
tmp1 = ao_integrals_n_e (i,j)
|
|
||||||
tmp2 = ao_integrals_n_e_cgtos(i,j)
|
|
||||||
dif = dabs(tmp1 - tmp2)
|
|
||||||
if(dif .gt. 1d-10) then
|
|
||||||
print*, ' error on:', i, j
|
|
||||||
print*, tmp1, tmp2, dif
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
acc += dif
|
|
||||||
nrm += dabs(tmp1)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
print *, ' acc (%) = ', 100.d0 * acc / nrm
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
|
|
||||||
subroutine check_ao_two_e_integral_cgtos()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i, j, k, l
|
|
||||||
double precision :: acc, nrm, dif
|
|
||||||
double precision :: tmp1, tmp2
|
|
||||||
double precision :: t1, t2, tt
|
|
||||||
|
|
||||||
double precision, external :: ao_two_e_integral
|
|
||||||
double precision, external :: ao_two_e_integral_cgtos
|
|
||||||
|
|
||||||
acc = 0.d0
|
|
||||||
nrm = 0.d0
|
|
||||||
|
|
||||||
tt = 0.d0
|
|
||||||
do i = 1, ao_num
|
|
||||||
!do i = 1, 1
|
|
||||||
call wall_time(t1)
|
|
||||||
do j = 1, ao_num
|
|
||||||
!do j = 1, 1
|
|
||||||
do k = 1, ao_num
|
|
||||||
!do k = 1, 1
|
|
||||||
do l = 1, ao_num
|
|
||||||
!do l = 21, 21
|
|
||||||
|
|
||||||
!call deb_ao_2eint_cgtos(i, j, k, l)
|
|
||||||
|
|
||||||
tmp1 = ao_two_e_integral (i, j, k, l)
|
|
||||||
tmp2 = ao_two_e_integral_cgtos(i, j, k, l)
|
|
||||||
|
|
||||||
dif = abs(tmp1 - tmp2)
|
|
||||||
if(dif .gt. 1d-10) then
|
|
||||||
print*, ' error on:', i, j, k, l
|
|
||||||
print*, tmp1, tmp2, dif
|
|
||||||
!stop
|
|
||||||
endif
|
|
||||||
acc += dif
|
|
||||||
nrm += abs(tmp1)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
call wall_time(t2)
|
|
||||||
tt += t2 - t1
|
|
||||||
print*, " % done = ", 100.d0 * dble(i) / ao_num
|
|
||||||
print*, ' ellapsed time (sec) =', tt
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!print *, ' acc (%) = ', 100.d0 * acc / nrm
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint1()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, n, i_rho
|
|
||||||
double precision :: dif_thr
|
|
||||||
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
|
|
||||||
complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), &
|
|
||||||
(+1d-9, +1d-6), &
|
|
||||||
(-1d-6, -1d-5), &
|
|
||||||
(+1d-3, -1d-2), &
|
|
||||||
(-1d-1, +1d-1), &
|
|
||||||
(+1d-0, +1d-1), &
|
|
||||||
(-1d+1, +1d+1), &
|
|
||||||
(+1d+2, +1d+1), &
|
|
||||||
(-1d+3, +1d+2), &
|
|
||||||
(+1d+4, +1d+4) /)
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int_an, int_nm
|
|
||||||
|
|
||||||
double precision, external :: rint
|
|
||||||
complex*16, external :: crint_1, crint_2
|
|
||||||
|
|
||||||
n = 10
|
|
||||||
dif_thr = 1d-7
|
|
||||||
|
|
||||||
do i_rho = 8, 10
|
|
||||||
!do i_rho = 7, 7
|
|
||||||
|
|
||||||
!rho = (-10.d0, 0.1d0)
|
|
||||||
!rho = (+10.d0, 0.1d0)
|
|
||||||
rho = rho_test(i_rho)
|
|
||||||
print*, "rho = ", real(rho), aimag(rho)
|
|
||||||
|
|
||||||
acc_re = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
acc_im = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i = 0, n
|
|
||||||
!int_an = crint_1(i, rho)
|
|
||||||
int_an = crint_2(i, rho)
|
|
||||||
call crint_quad_1(i, rho, 100000000, int_nm)
|
|
||||||
|
|
||||||
dif_re = dabs(real(int_an) - real(int_nm))
|
|
||||||
dif_im = dabs(aimag(int_an) - aimag(int_nm))
|
|
||||||
|
|
||||||
if((dif_re .gt. dif_thr) .or. (dif_im .gt. dif_thr)) then
|
|
||||||
print*, ' error on i =', i
|
|
||||||
print*, real(int_an), real(int_nm), dif_re
|
|
||||||
print*, aimag(int_an), aimag(int_nm), dif_im
|
|
||||||
!print*, rint(i, real(rho))
|
|
||||||
print*, crint_1(i, rho)
|
|
||||||
!print*, crint_2(i, rho)
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
acc_re += dif_re
|
|
||||||
nrm_re += dabs(real(int_nm))
|
|
||||||
acc_im += dif_im
|
|
||||||
nrm_im += dabs(aimag(int_nm))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15)
|
|
||||||
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint2()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i, n, i_rho
|
|
||||||
double precision :: dif_thr
|
|
||||||
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
|
|
||||||
complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), &
|
|
||||||
(+1d-9, +1d-6), &
|
|
||||||
(-1d-6, -1d-5), &
|
|
||||||
(+1d-3, -1d-2), &
|
|
||||||
(-1d-1, +1d-1), &
|
|
||||||
(+1d-0, +1d-1), &
|
|
||||||
(-1d+1, +1d+1), &
|
|
||||||
(+1d+2, +1d+1), &
|
|
||||||
(-1d+3, +1d+2), &
|
|
||||||
(+1d+4, +1d+4) /)
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int_an, int_nm
|
|
||||||
complex*16, external :: crint_1, crint_2
|
|
||||||
|
|
||||||
n = 30
|
|
||||||
dif_thr = 1d-12
|
|
||||||
|
|
||||||
do i_rho = 1, 10
|
|
||||||
rho = rho_test(i_rho)
|
|
||||||
print*, "rho = ", real(rho), aimag(rho)
|
|
||||||
|
|
||||||
acc_re = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
acc_im = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i = 0, n
|
|
||||||
int_an = crint_1(i, rho)
|
|
||||||
int_nm = crint_2(i, rho)
|
|
||||||
|
|
||||||
dif_re = dabs(real(int_an) - real(int_nm))
|
|
||||||
!if(dif_re .gt. dif_thr) then
|
|
||||||
! print*, ' error in real part:', i
|
|
||||||
! print*, real(int_an), real(int_nm), dif_re
|
|
||||||
! stop
|
|
||||||
!endif
|
|
||||||
acc_re += dif_re
|
|
||||||
nrm_re += dabs(real(int_nm))
|
|
||||||
|
|
||||||
dif_im = dabs(aimag(int_an) - aimag(int_nm))
|
|
||||||
!if(dif_im .gt. dif_thr) then
|
|
||||||
! print*, ' error in imag part:', i
|
|
||||||
! print*, aimag(int_an), aimag(int_nm), dif_im
|
|
||||||
! stop
|
|
||||||
!endif
|
|
||||||
acc_im += dif_im
|
|
||||||
nrm_im += dabs(aimag(int_nm))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15)
|
|
||||||
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint3()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i_test, n_test
|
|
||||||
integer :: nx, ny, n, n_quad
|
|
||||||
integer :: i, seed_size, clock_time
|
|
||||||
double precision :: xr(1:4), x
|
|
||||||
double precision :: yr(1:4), y
|
|
||||||
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
|
|
||||||
double precision :: delta_ref
|
|
||||||
double precision :: t1, t2, t_int1, t_int2
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int1_old, int1_ref, int2_old, int2_ref
|
|
||||||
integer, allocatable :: seed(:)
|
|
||||||
|
|
||||||
complex*16, external :: crint_2
|
|
||||||
|
|
||||||
call random_seed(size=seed_size)
|
|
||||||
allocate(seed(seed_size))
|
|
||||||
call system_clock(count=clock_time)
|
|
||||||
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
|
|
||||||
!seed = 123456789
|
|
||||||
call random_seed(put=seed)
|
|
||||||
|
|
||||||
|
|
||||||
t_int1 = 0.d0
|
|
||||||
t_int2 = 0.d0
|
|
||||||
|
|
||||||
n_test = 1
|
|
||||||
|
|
||||||
acc_re = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
acc_im = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i_test = 1, n_test
|
|
||||||
|
|
||||||
! Re(rho)
|
|
||||||
call random_number(xr)
|
|
||||||
x = xr(1)
|
|
||||||
if(xr(2) .gt. 0.5d0) x = -x
|
|
||||||
nx = int(15.d0 * xr(3))
|
|
||||||
if(xr(4) .gt. 0.5d0) nx = -nx
|
|
||||||
x = x * 10.d0**nx
|
|
||||||
|
|
||||||
! Im(rho)
|
|
||||||
call random_number(yr)
|
|
||||||
y = yr(1)
|
|
||||||
if(yr(2) .gt. 0.5d0) y = -y
|
|
||||||
ny = int(5.d0 * yr(3))
|
|
||||||
if(yr(4) .gt. 0.5d0) ny = -ny
|
|
||||||
y = y * 10.d0**ny
|
|
||||||
|
|
||||||
rho = x + (0.d0, 1.d0) * y
|
|
||||||
|
|
||||||
call random_number(x)
|
|
||||||
x = 31.d0 * x
|
|
||||||
n = int(x)
|
|
||||||
!if(n.eq.0) cycle
|
|
||||||
|
|
||||||
n = 0
|
|
||||||
!rho = (-6.83897018210218d0, -7.24479852507338d0)
|
|
||||||
rho = (-9.83206247355480d0, 0.445269582329036d0)
|
|
||||||
|
|
||||||
print*, " n = ", n
|
|
||||||
print*, " rho = ", real(rho), aimag(rho)
|
|
||||||
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
int1_old = crint_2(n, rho)
|
|
||||||
!n_quad = 10000000
|
|
||||||
!call crint_quad_1(n, rho, n_quad, int1_old)
|
|
||||||
!!delta_ref = 1.d0
|
|
||||||
!!do while(delta_ref .gt. 1d-12)
|
|
||||||
!! n_quad = n_quad * 10
|
|
||||||
!! !print*, " delta = ", delta_ref
|
|
||||||
!! !print*, " increasing n_quad to:", n_quad
|
|
||||||
!! call crint_quad_1(n, rho, n_quad, int1_ref)
|
|
||||||
!! delta_ref = abs(int1_ref - int1_old)
|
|
||||||
!! int1_old = int1_ref
|
|
||||||
!! if(n_quad .ge. 1000000000) then
|
|
||||||
!! print*, ' convergence was not reached for crint_quad_1'
|
|
||||||
!! print*, " delta = ", delta_ref
|
|
||||||
!! exit
|
|
||||||
!! endif
|
|
||||||
!!enddo
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int1 = t_int1 + t2 - t1
|
|
||||||
!print*, " n_quad for crint_quad_1:", n_quad
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
n_quad = 10000000
|
|
||||||
call crint_quad_12(n, rho, n_quad, int2_old)
|
|
||||||
!delta_ref = 1.d0
|
|
||||||
!do while(delta_ref .gt. 1d-12)
|
|
||||||
! n_quad = n_quad * 10
|
|
||||||
! !print*, " delta = ", delta_ref
|
|
||||||
! !print*, " increasing n_quad to:", n_quad
|
|
||||||
! call crint_quad_12(n, rho, n_quad, int2_ref)
|
|
||||||
! delta_ref = abs(int2_ref - int2_old)
|
|
||||||
! int2_old = int2_ref
|
|
||||||
! if(n_quad .ge. 1000000000) then
|
|
||||||
! print*, ' convergence was not reached for crint_quad_2'
|
|
||||||
! print*, " delta = ", delta_ref
|
|
||||||
! exit
|
|
||||||
! endif
|
|
||||||
!enddo
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int2 = t_int2 + t2 - t1
|
|
||||||
!print*, " n_quad for crint_quad_2:", n_quad
|
|
||||||
|
|
||||||
dif_re = dabs(real(int1_old) - real(int2_old))
|
|
||||||
dif_im = dabs(aimag(int1_old) - aimag(int2_old))
|
|
||||||
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
|
|
||||||
print*, ' important error found: '
|
|
||||||
print*, " n = ", n
|
|
||||||
print*, " rho = ", real(rho), aimag(rho)
|
|
||||||
print*, real(int1_old), real(int2_old), dif_re
|
|
||||||
print*, aimag(int1_old), aimag(int2_old), dif_im
|
|
||||||
!stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
if((real(int1_old) /= real(int1_old)) .or. (aimag(int1_old) /= aimag(int1_old)) .or. &
|
|
||||||
(real(int2_old) /= real(int2_old)) .or. (aimag(int2_old) /= aimag(int2_old)) ) then
|
|
||||||
cycle
|
|
||||||
else
|
|
||||||
acc_re += dif_re
|
|
||||||
acc_im += dif_im
|
|
||||||
nrm_re += dabs(real(int1_old))
|
|
||||||
nrm_im += dabs(aimag(int1_old))
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15)
|
|
||||||
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15)
|
|
||||||
|
|
||||||
print*, "crint_quad_1 wall time (sec) = ", t_int1
|
|
||||||
print*, "crint_quad_2 wall time (sec) = ", t_int2
|
|
||||||
|
|
||||||
|
|
||||||
deallocate(seed)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint4()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i_test, n_test
|
|
||||||
integer :: i, seed_size, clock_time
|
|
||||||
double precision :: xr(1), x, shift
|
|
||||||
double precision :: yr(1), y
|
|
||||||
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
|
|
||||||
double precision :: t1, t2, t_int1, t_int2
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int1, int2, int3
|
|
||||||
integer, allocatable :: seed(:)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
call random_seed(size=seed_size)
|
|
||||||
allocate(seed(seed_size))
|
|
||||||
call system_clock(count=clock_time)
|
|
||||||
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
|
|
||||||
!seed = 123456789
|
|
||||||
call random_seed(put=seed)
|
|
||||||
|
|
||||||
|
|
||||||
t_int1 = 0.d0
|
|
||||||
t_int2 = 0.d0
|
|
||||||
|
|
||||||
n_test = 100
|
|
||||||
shift = 15.d0
|
|
||||||
|
|
||||||
acc_re = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
acc_im = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i_test = 1, n_test
|
|
||||||
|
|
||||||
call random_number(xr)
|
|
||||||
call random_number(yr)
|
|
||||||
|
|
||||||
x = 1.d0 * (2.d0 * shift * xr(1) - shift)
|
|
||||||
y = 1.d0 * (2.d0 * shift * yr(1) - shift)
|
|
||||||
|
|
||||||
rho = x + (0.d0, 1.d0) * y
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
call zboysfun00_1(rho, int1)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int1 = t_int1 + t2 - t1
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
call zboysfun00_2(rho, int2)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int2 = t_int2 + t2 - t1
|
|
||||||
|
|
||||||
dif_re = dabs(real(int1) - real(int2))
|
|
||||||
dif_im = dabs(aimag(int1) - aimag(int2))
|
|
||||||
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
|
|
||||||
print*, ' important error found: '
|
|
||||||
print*, " rho = ", x, y
|
|
||||||
print*, real(int1), real(int2), dif_re
|
|
||||||
print*, aimag(int1), aimag(int2), dif_im
|
|
||||||
call crint_quad_12(0, rho, 10000000, int3)
|
|
||||||
if(zabs(int1 - int3) .lt. zabs(int2 - int3)) then
|
|
||||||
print*, ' implementation 2 seems to be wrong'
|
|
||||||
else
|
|
||||||
print*, ' implementation 1 seems to be wrong'
|
|
||||||
print*, ' quad 10000000:', real(int3), aimag(int3)
|
|
||||||
call crint_quad_12(0, rho, 100000000, int3)
|
|
||||||
print*, ' quad 100000000:', real(int3), aimag(int3)
|
|
||||||
endif
|
|
||||||
!print*, ' quad:', real(int3), aimag(int3)
|
|
||||||
!stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
if((real(int1) /= real(int1)) .or. (aimag(int1) /= aimag(int1)) .or. &
|
|
||||||
(real(int2) /= real(int2)) .or. (aimag(int2) /= aimag(int2)) ) then
|
|
||||||
cycle
|
|
||||||
else
|
|
||||||
acc_re += dif_re
|
|
||||||
acc_im += dif_im
|
|
||||||
nrm_re += dabs(real(int1))
|
|
||||||
nrm_im += dabs(aimag(int1))
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15)
|
|
||||||
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15)
|
|
||||||
|
|
||||||
print*, "zerf_1 wall time (sec) = ", t_int1
|
|
||||||
print*, "zerf_2 wall time (sec) = ", t_int2
|
|
||||||
|
|
||||||
|
|
||||||
deallocate(seed)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint5()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i_test, n_test
|
|
||||||
integer :: i, seed_size, clock_time
|
|
||||||
integer :: n
|
|
||||||
double precision :: xr(1), yr(1), nr(1), x, shift, y
|
|
||||||
double precision :: dif1_re, dif1_im, acc1_re, acc1_im
|
|
||||||
double precision :: dif2_re, dif2_im, acc2_re, acc2_im
|
|
||||||
double precision :: nrm_re, nrm_im
|
|
||||||
double precision :: t1, t2, t_int1, t_int2
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int1, int2, int_ref
|
|
||||||
integer, allocatable :: seed(:)
|
|
||||||
|
|
||||||
complex*16, external :: crint_1, crint_2
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
call random_seed(size=seed_size)
|
|
||||||
allocate(seed(seed_size))
|
|
||||||
call system_clock(count=clock_time)
|
|
||||||
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
|
|
||||||
!seed = 123456789
|
|
||||||
call random_seed(put=seed)
|
|
||||||
|
|
||||||
|
|
||||||
t_int1 = 0.d0
|
|
||||||
t_int2 = 0.d0
|
|
||||||
|
|
||||||
n_test = 100
|
|
||||||
|
|
||||||
acc1_re = 0.d0
|
|
||||||
acc1_im = 0.d0
|
|
||||||
acc2_re = 0.d0
|
|
||||||
acc2_im = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i_test = 1, n_test
|
|
||||||
|
|
||||||
call random_number(xr)
|
|
||||||
call random_number(yr)
|
|
||||||
call random_number(nr)
|
|
||||||
|
|
||||||
x = 1.d+1 * (30.d0 * xr(1) - 15.d0)
|
|
||||||
y = 1.d+1 * (30.d0 * yr(1) - 15.d0)
|
|
||||||
n = int(16.d0 * nr(1))
|
|
||||||
|
|
||||||
rho = x + (0.d0, 1.d0) * y
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
int1 = crint_1(n, rho)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int1 = t_int1 + t2 - t1
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
int2 = crint_2(n, rho)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int2 = t_int2 + t2 - t1
|
|
||||||
|
|
||||||
call crint_quad_12(n, rho, 10000000, int_ref)
|
|
||||||
|
|
||||||
dif1_re = dabs(real(int1) - real(int_ref))
|
|
||||||
dif1_im = dabs(aimag(int1) - aimag(int_ref))
|
|
||||||
|
|
||||||
dif2_re = dabs(real(int2) - real(int_ref))
|
|
||||||
dif2_im = dabs(aimag(int2) - aimag(int_ref))
|
|
||||||
|
|
||||||
if((dif2_re .gt. 1d-7) .or. (dif2_im .gt. 1d-7)) then
|
|
||||||
print*, ' important error found: '
|
|
||||||
print*, " n, rho = ", n, x, y
|
|
||||||
print*, real(int1), real(int2), real(int_ref)
|
|
||||||
print*, aimag(int1), aimag(int2), aimag(int_ref)
|
|
||||||
!stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
acc1_re += dif1_re
|
|
||||||
acc1_im += dif1_im
|
|
||||||
|
|
||||||
acc2_re += dif2_re
|
|
||||||
acc2_im += dif2_im
|
|
||||||
|
|
||||||
nrm_re += dabs(real(int_ref))
|
|
||||||
nrm_im += dabs(aimag(int_ref))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "accuracy on boys_1 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc1_im / (nrm_im + 1d-15)
|
|
||||||
print*, "accuracy on boys_2 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc2_im / (nrm_im + 1d-15)
|
|
||||||
|
|
||||||
print*, "boys_1 wall time (sec) = ", t_int1
|
|
||||||
print*, "boys_2 wall time (sec) = ", t_int2
|
|
||||||
|
|
||||||
|
|
||||||
deallocate(seed)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine check_crint6()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer :: i_test, n_test
|
|
||||||
integer :: i, seed_size, clock_time
|
|
||||||
integer :: n
|
|
||||||
double precision :: xr(1), yr(1), nr(1), x, shift, y
|
|
||||||
double precision :: dif_re, dif_im, acc_re, acc_im
|
|
||||||
double precision :: nrm_re, nrm_im
|
|
||||||
double precision :: t1, t2, t_int1, t_int2
|
|
||||||
complex*16 :: rho
|
|
||||||
complex*16 :: int1, int2, int3
|
|
||||||
integer, allocatable :: seed(:)
|
|
||||||
|
|
||||||
complex*16, external :: crint_1, crint_2
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
call random_seed(size=seed_size)
|
|
||||||
allocate(seed(seed_size))
|
|
||||||
call system_clock(count=clock_time)
|
|
||||||
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
|
|
||||||
!seed = 123456789
|
|
||||||
call random_seed(put=seed)
|
|
||||||
|
|
||||||
|
|
||||||
t_int1 = 0.d0
|
|
||||||
t_int2 = 0.d0
|
|
||||||
|
|
||||||
n_test = 100
|
|
||||||
|
|
||||||
acc_re = 0.d0
|
|
||||||
acc_im = 0.d0
|
|
||||||
nrm_re = 0.d0
|
|
||||||
nrm_im = 0.d0
|
|
||||||
do i_test = 1, n_test
|
|
||||||
|
|
||||||
call random_number(xr)
|
|
||||||
call random_number(yr)
|
|
||||||
call random_number(nr)
|
|
||||||
|
|
||||||
x = 1.d0 * (30.d0 * xr(1) - 15.d0)
|
|
||||||
y = 1.d0 * (30.d0 * yr(1) - 15.d0)
|
|
||||||
n = int(16.d0 * nr(1))
|
|
||||||
|
|
||||||
rho = x + (0.d0, 1.d0) * y
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
int1 = crint_1(n, rho)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int1 = t_int1 + t2 - t1
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
int2 = crint_2(n, rho)
|
|
||||||
call wall_time(t2)
|
|
||||||
t_int2 = t_int2 + t2 - t1
|
|
||||||
|
|
||||||
dif_re = dabs(real(int1) - real(int2))
|
|
||||||
dif_im = dabs(aimag(int1) - aimag(int2))
|
|
||||||
|
|
||||||
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
|
|
||||||
print*, ' important error found: '
|
|
||||||
print*, " n, rho = ", n, x, y
|
|
||||||
print*, real(int1), real(int2), dif_re
|
|
||||||
print*, aimag(int1), aimag(int2), dif_im
|
|
||||||
call crint_quad_12(n, rho, 100000000, int3)
|
|
||||||
print*, ' quad 100000000:', real(int3), aimag(int3)
|
|
||||||
!print*, ' quad 100000000:', dabs(real(int1) - real(int3)), dabs(aimag(int1) - aimag(int3))
|
|
||||||
!stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
acc_re += dif_re
|
|
||||||
acc_im += dif_im
|
|
||||||
nrm_re += dabs(real(int1))
|
|
||||||
nrm_im += dabs(aimag(int1))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print*, "diff (%):", 100.d0 * acc_re / (nrm_re + 1d-15), 100.d0 * acc_im / (nrm_im + 1d-15)
|
|
||||||
|
|
||||||
print*, "boys_1 wall time (sec) = ", t_int1
|
|
||||||
print*, "boys_2 wall time (sec) = ", t_int2
|
|
||||||
|
|
||||||
|
|
||||||
deallocate(seed)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine deb_ao_2eint_cgtos(i, j, k, l)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! integral of the AO basis <ik|jl> or (ij|kl)
|
|
||||||
! i(r1) j(r1) 1/r12 k(r2) l(r2)
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'utils/constants.include.F'
|
|
||||||
|
|
||||||
integer, intent(in) :: i, j, k, l
|
|
||||||
|
|
||||||
integer :: p, q, r, s
|
|
||||||
integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3)
|
|
||||||
integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3)
|
|
||||||
complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3)
|
|
||||||
complex*16 :: expo1, expo2, expo3, expo4
|
|
||||||
complex*16 :: P1_center(3), pp1
|
|
||||||
complex*16 :: P2_center(3), pp2
|
|
||||||
complex*16 :: Q1_center(3), qq1
|
|
||||||
complex*16 :: Q2_center(3), qq2
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
dim1 = n_pt_max_integrals
|
|
||||||
|
|
||||||
num_i = ao_nucl(i)
|
|
||||||
num_j = ao_nucl(j)
|
|
||||||
num_k = ao_nucl(k)
|
|
||||||
num_l = ao_nucl(l)
|
|
||||||
|
|
||||||
if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then
|
|
||||||
|
|
||||||
!print*, ao_prim_num(i), ao_prim_num(j), ao_prim_num(k), ao_prim_num(l)
|
|
||||||
|
|
||||||
do p = 1, 3
|
|
||||||
I_power(p) = ao_power(i,p)
|
|
||||||
J_power(p) = ao_power(j,p)
|
|
||||||
K_power(p) = ao_power(k,p)
|
|
||||||
L_power(p) = ao_power(l,p)
|
|
||||||
I_center(p) = nucl_coord(num_i,p) * (1.d0, 0.d0)
|
|
||||||
J_center(p) = nucl_coord(num_j,p) * (1.d0, 0.d0)
|
|
||||||
K_center(p) = nucl_coord(num_k,p) * (1.d0, 0.d0)
|
|
||||||
L_center(p) = nucl_coord(num_l,p) * (1.d0, 0.d0)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do p = 1, ao_prim_num(i)
|
|
||||||
expo1 = ao_expo_cgtos_ord_transp(p,i)
|
|
||||||
!print*, "expo1 = ", expo1
|
|
||||||
!print*, "center1 = ", I_center
|
|
||||||
|
|
||||||
do q = 1, ao_prim_num(j)
|
|
||||||
expo2 = ao_expo_cgtos_ord_transp(q,j)
|
|
||||||
!print*, "expo2 = ", expo2
|
|
||||||
!print*, "center2 = ", J_center
|
|
||||||
|
|
||||||
pp1 = expo1 + expo2
|
|
||||||
P1_center(1:3) = (expo1 * I_center(1:3) + expo2 * J_center(1:3)) / pp1
|
|
||||||
iorder_p1(1:3) = I_power(1:3) + J_power(1:3)
|
|
||||||
|
|
||||||
pp2 = conjg(expo1) + expo2
|
|
||||||
P2_center(1:3) = (conjg(expo1) * I_center(1:3) + expo2 * J_center(1:3)) / pp2
|
|
||||||
iorder_p2(1:3) = I_power(1:3) + J_power(1:3)
|
|
||||||
|
|
||||||
do r = 1, ao_prim_num(k)
|
|
||||||
expo3 = ao_expo_cgtos_ord_transp(r,k)
|
|
||||||
!print*, "expo3 = ", expo3
|
|
||||||
!print*, "center3 = ", K_center
|
|
||||||
|
|
||||||
do s = 1, ao_prim_num(l)
|
|
||||||
expo4 = ao_expo_cgtos_ord_transp(s,l)
|
|
||||||
!print*, "expo4 = ", expo4
|
|
||||||
!print*, "center4 = ", L_center
|
|
||||||
|
|
||||||
qq1 = expo3 + expo4
|
|
||||||
Q1_center(1:3) = (expo3 * K_center(1:3) + expo4 * L_center(1:3)) / qq1
|
|
||||||
iorder_q1(1:3) = K_power(1:3) + L_power(1:3)
|
|
||||||
|
|
||||||
qq2 = conjg(expo3) + expo4
|
|
||||||
Q2_center(1:3) = (conjg(expo3) * K_center(1:3) + expo4 * L_center(1:3)) / qq2
|
|
||||||
iorder_q2(1:3) = K_power(1:3) + L_power(1:3)
|
|
||||||
|
|
||||||
call deb_cboys(P1_center, pp1, iorder_p1, Q1_center, qq1, iorder_q1)
|
|
||||||
call deb_cboys(P1_center, pp1, iorder_p1, Q2_center, qq2, iorder_q2)
|
|
||||||
call deb_cboys(P2_center, pp2, iorder_p2, Q1_center, qq1, iorder_q1)
|
|
||||||
call deb_cboys(P2_center, pp2, iorder_p2, Q2_center, qq2, iorder_q2)
|
|
||||||
call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q1_center, qq1, iorder_q1)
|
|
||||||
call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q2_center, qq2, iorder_q2)
|
|
||||||
call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q1_center, qq1, iorder_q1)
|
|
||||||
call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q2_center, qq2, iorder_q2)
|
|
||||||
enddo ! s
|
|
||||||
enddo ! r
|
|
||||||
enddo ! q
|
|
||||||
enddo ! p
|
|
||||||
|
|
||||||
endif ! same centers
|
|
||||||
|
|
||||||
return
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine deb_cboys(P_center, p, iorder_p, Q_center, q, iorder_q)
|
|
||||||
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'utils/constants.include.F'
|
|
||||||
|
|
||||||
integer, intent(in) :: iorder_p(3), iorder_q(3)
|
|
||||||
complex*16, intent(in) :: P_center(3), p
|
|
||||||
complex*16, intent(in) :: Q_center(3), q
|
|
||||||
|
|
||||||
integer :: iorder, n
|
|
||||||
complex*16 :: dist, rho
|
|
||||||
complex*16 :: int1, int2
|
|
||||||
|
|
||||||
complex*16, external :: crint_2
|
|
||||||
|
|
||||||
|
|
||||||
dist = (P_center(1) - Q_center(1)) * (P_center(1) - Q_center(1)) &
|
|
||||||
+ (P_center(2) - Q_center(2)) * (P_center(2) - Q_center(2)) &
|
|
||||||
+ (P_center(3) - Q_center(3)) * (P_center(3) - Q_center(3))
|
|
||||||
rho = dist * p * q / (p + q)
|
|
||||||
|
|
||||||
if(real(rho) .lt. -5.d0) then
|
|
||||||
print*, 'warning ! impotant negative rho: ', rho
|
|
||||||
endif
|
|
||||||
|
|
||||||
!if(abs(rho) .lt. 1d-15) return
|
|
||||||
|
|
||||||
iorder = 2*iorder_p(1)+2*iorder_q(1) + 2*iorder_p(2)+2*iorder_q(2) + 2*iorder_p(3)+2*iorder_q(3)
|
|
||||||
n = shiftr(iorder, 1)
|
|
||||||
|
|
||||||
!write(33,*) n, real(rho), aimag(rho)
|
|
||||||
!print*, n, real(rho), aimag(rho)
|
|
||||||
|
|
||||||
int1 = crint_2(n, rho)
|
|
||||||
call crint_quad_12(n, rho, 1000000, int2)
|
|
||||||
|
|
||||||
if(abs(int1 - int2) .gt. 1d-5) then
|
|
||||||
print*, ' important error found: '
|
|
||||||
print*, p!, P_center
|
|
||||||
print*, q!, Q_center
|
|
||||||
print*, dist
|
|
||||||
print*, " n, tho = ", n, real(rho), aimag(rho)
|
|
||||||
print*, real(int1), real(int2), dabs(real(int1-int2))
|
|
||||||
print*, aimag(int1), aimag(int2), dabs(aimag(int1-int2))
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue
Block a user