9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-09 15:03:37 +01:00
qp2/plugins/local/non_h_ints_mu/debug_fit.irp.f

450 lines
11 KiB
Fortran
Raw Normal View History

2023-02-06 19:00:35 +01:00
! --
program debug_fit
implicit none
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-02-06 19:00:35 +01:00
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
2024-01-15 12:02:38 +01:00
PROVIDE j2e_type mu_erf
PROVIDE j1e_type j1e_coef j1e_expo
PROVIDE env_type env_coef env_expo
provide tc_integ_type
2023-02-06 19:00:35 +01:00
2024-01-15 12:02:38 +01:00
if(tc_integ_type .eq. "numeric") then
2023-09-05 11:52:08 +02:00
my_extra_grid_becke = .True.
PROVIDE tc_grid2_a tc_grid2_r
my_n_pt_r_extra_grid = tc_grid2_r
my_n_pt_a_extra_grid = tc_grid2_a
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
endif
2024-01-15 12:02:38 +01:00
!call test_env_nucl()
!call test_grad_env_nucl()
2023-02-06 19:00:35 +01:00
2023-05-08 23:31:20 +02:00
!call test_fit_u()
2023-02-06 19:00:35 +01:00
!call test_fit_u2()
!call test_fit_ugradu()
2023-09-05 11:52:08 +02:00
call test_grad1_u12_withsq_num()
2023-02-06 19:00:35 +01:00
end
! ---
2024-01-15 12:02:38 +01:00
subroutine test_env_nucl()
2023-02-06 19:00:35 +01:00
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
2024-01-15 12:02:38 +01:00
double precision, external :: env_nucl
2023-02-06 19:00:35 +01:00
2024-01-15 12:02:38 +01:00
print*, ' test_env_nucl ...'
2023-02-06 19:00:35 +01:00
2024-01-15 12:02:38 +01:00
PROVIDE env_val
2023-02-06 19:00:35 +01:00
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
2024-01-15 12:02:38 +01:00
i_exc = env_val(ipoint)
i_num = env_nucl(r)
2023-02-06 19:00:35 +01:00
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
2024-01-15 12:02:38 +01:00
print *, ' problem in env_val on', ipoint
2023-02-06 19:00:35 +01:00
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
2024-01-15 12:02:38 +01:00
end
2023-02-06 19:00:35 +01:00
! ---
2024-01-15 12:02:38 +01:00
subroutine test_grad_env_nucl()
2023-02-06 19:00:35 +01:00
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
2024-01-15 12:02:38 +01:00
double precision, external :: grad_x_env_nucl_num
double precision, external :: grad_y_env_nucl_num
double precision, external :: grad_z_env_nucl_num
2023-02-06 19:00:35 +01:00
2024-01-15 12:02:38 +01:00
PROVIDE env_grad
2023-02-06 19:00:35 +01:00
2024-01-15 12:02:38 +01:00
print*, ' test_grad_env_nucl ...'
2023-02-06 19:00:35 +01:00
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
2024-01-15 12:02:38 +01:00
i_exc = env_grad(1,ipoint)
i_num = grad_x_env_nucl_num(r)
2023-02-06 19:00:35 +01:00
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
2024-01-15 12:02:38 +01:00
print *, ' problem in x of env_grad on', ipoint
2023-02-06 19:00:35 +01:00
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
2024-01-15 12:02:38 +01:00
i_exc = env_grad(2,ipoint)
i_num = grad_y_env_nucl_num(r)
2023-02-06 19:00:35 +01:00
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
2024-01-15 12:02:38 +01:00
print *, ' problem in y of env_grad on', ipoint
2023-02-06 19:00:35 +01:00
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
2024-01-15 12:02:38 +01:00
i_exc = env_grad(3,ipoint)
i_num = grad_z_env_nucl_num(r)
2023-02-06 19:00:35 +01:00
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
2024-01-15 12:02:38 +01:00
print *, ' problem in z of env_grad on', ipoint
2023-02-06 19:00:35 +01:00
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
2024-01-15 12:02:38 +01:00
end
2023-02-06 19:00:35 +01:00
! ---
subroutine test_fit_ugradu()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2, tmp, dx, dy, dz
double precision :: r1(3), r2(3), grad(3)
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_ugradu ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_1_erf(i)
coef = coef_gauss_j_mu_1_erf(i)
i_fit += coef * dexp(-expo*x2)
enddo
i_fit = i_fit / dsqrt(x2)
tmp = j12_mu(r1, r2)
2023-04-21 20:50:00 +02:00
call grad1_j12_mu(r1, r2, grad)
2023-02-06 19:00:35 +01:00
! ---
i_exc = tmp * grad(1)
i_num = i_fit * dx
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on x in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = tmp * grad(2)
i_num = i_fit * dy
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on y in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = tmp * grad(3)
i_num = i_fit * dz
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on z in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
2024-01-15 12:02:38 +01:00
end
2023-02-06 19:00:35 +01:00
! ---
subroutine test_fit_u()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3), dx, dy, dz
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_u ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x(i)
coef = coef_gauss_j_mu_x(i)
i_fit += coef * dexp(-expo*x2)
enddo
i_exc = j12_mu(r1, r2)
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
2024-01-15 12:02:38 +01:00
end
2023-02-06 19:00:35 +01:00
! ---
subroutine test_fit_u2()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3), dx, dy, dz, tmp
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_u2 ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x_2(i)
coef = coef_gauss_j_mu_x_2(i)
i_fit += coef * dexp(-expo*x2)
enddo
tmp = j12_mu(r1, r2)
i_exc = tmp * tmp
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
2024-01-15 12:02:38 +01:00
end
2023-02-06 19:00:35 +01:00
! ---
2023-09-05 11:52:08 +02:00
subroutine test_grad1_u12_withsq_num()
implicit none
integer :: ipoint, jpoint, m
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision, allocatable :: tmp_grad1_u12_squared(:,:), tmp_grad1_u12(:,:,:)
print*, ' test_grad1_u12_withsq_num ...'
PROVIDE grad1_u12_num grad1_u12_squared_num
allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_points_final_grid))
allocate(tmp_grad1_u12(n_points_extra_final_grid,n_points_final_grid,3))
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
2024-01-23 13:25:16 +01:00
call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, tmp_grad1_u12(1,ipoint,1) &
, tmp_grad1_u12(1,ipoint,2) &
, tmp_grad1_u12(1,ipoint,3) &
, tmp_grad1_u12_squared(1,ipoint))
2023-09-05 11:52:08 +02:00
do jpoint = 1, n_points_extra_final_grid
i_exc = grad1_u12_squared_num(jpoint,ipoint)
i_num = tmp_grad1_u12_squared(jpoint,ipoint)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in grad1_u12_squared_num on', ipoint, jpoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
stop
endif
acc_tot += acc_ij
normalz += dabs(i_num)
do m = 1, 3
i_exc = grad1_u12_num(jpoint,ipoint,m)
i_num = tmp_grad1_u12(jpoint,ipoint,m)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in grad1_u12_num on', ipoint, jpoint, m
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
stop
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
enddo
!print*, ' acc_tot = ', acc_tot
!print*, ' normalz = ', normalz
print*, ' accuracy (%) = ', 100.d0 * acc_tot / normalz
return
2024-01-15 12:02:38 +01:00
end
2023-09-05 11:52:08 +02:00
! ---
2023-02-06 19:00:35 +01:00