LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-01-03 17:15:40 +01:00
Code Releases Activity

it compiles !

Browse Source
This commit is contained in:
eginer 2022-10-22 19:09:22 +02:00
parent fdc5b7a467
commit aa55f7de28
1 changed files with 0 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
src/non_h_ints_mu/test_non_h_ints.irp.f
View File

@ -1,135 +0,0 @@
program test_non_h
implicit none
my_grid_becke = .True.
! my_n_pt_r_grid = 50
! my_n_pt_a_grid = 74
my_n_pt_r_grid = 10 ! small grid for quick debug
my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
!call routine_grad_squared
! call routine_fit
call routine_grad_squared_new
end
subroutine routine_lapl_grad
implicit none
integer :: i,j,k,l
double precision :: grad_lapl, get_ao_tc_sym_two_e_pot,new,accu,contrib
double precision :: ao_two_e_integral_erf,get_ao_two_e_integral,count_n,accu_relat
! !!!!!!!!!!!!!!!!!!!!! WARNING
! THIS ROUTINE MAKES SENSE ONLY IF HAND MODIFIED coef_gauss_eff_pot(1:n_max_fit_slat) = 0. to cancel (1-erf(mu*r12))^2
accu = 0.d0
accu_relat = 0.d0
count_n = 0.d0
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
grad_lapl = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) ! pure gaussian part : comes from Lapl
grad_lapl += ao_two_e_integral_erf(i, k, j, l) ! erf(mu r12)/r12 : comes from Lapl
grad_lapl += ao_non_hermit_term_chemist(k,i,l,j) ! \grad u(r12) . grad
new = tc_grad_and_lapl_ao(k,i,l,j)
new += get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
contrib = dabs(new - grad_lapl)
if(dabs(grad_lapl).gt.1.d-12)then
count_n += 1.d0
accu_relat += 2.0d0 * contrib/dabs(grad_lapl+new)
endif
if(contrib.gt.1.d-10)then
print*,i,j,k,l
print*,grad_lapl,new,contrib
print*,2.0d0*contrib/dabs(grad_lapl+new+1.d-12)
endif
accu += contrib
enddo
enddo
enddo
enddo
print*,'accu = ',accu/count_n
print*,'accu/rel = ',accu_relat/count_n
end
subroutine routine_grad_squared
implicit none
integer :: i,j,k,l
double precision :: grad_squared, get_ao_tc_sym_two_e_pot,new,accu,contrib
double precision :: count_n,accu_relat
! !!!!!!!!!!!!!!!!!!!!! WARNING
! THIS ROUTINE MAKES SENSE ONLY IF HAND MODIFIED coef_gauss_eff_pot(n_max_fit_slat:n_max_fit_slat+1) = 0. to cancel exp(-'mu*r12)^2)
accu = 0.d0
accu_relat = 0.d0
count_n = 0.d0
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
grad_squared = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) ! pure gaussian part : comes from Lapl
new = tc_grad_square_ao(k,i,l,j)
contrib = dabs(new - grad_squared)
if(dabs(grad_squared).gt.1.d-12)then
count_n += 1.d0
accu_relat += 2.0d0 * contrib/dabs(grad_squared+new)
endif
if(contrib.gt.1.d-10)then
print*,i,j,k,l
print*,grad_squared,new,contrib
print*,2.0d0*contrib/dabs(grad_squared+new+1.d-12)
endif
accu += contrib
enddo
enddo
enddo
enddo
print*,'accu = ',accu/count_n
print*,'accu/rel = ',accu_relat/count_n
end
subroutine routine_grad_squared_new
implicit none
integer :: i,j,k,l,ipoint
double precision :: grad_squared, get_ao_tc_sym_two_e_pot,new,accu,contrib
double precision :: count_n,accu_relat
accu = 0.d0
accu_relat = 0.d0
count_n = 0.d0
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
grad_squared = grad_1_squared_u_ij_mu(j,i,ipoint)
new = grad_1_squared_u_ij_mu_new(ipoint,j,i)
contrib = dabs(new - grad_squared)
if(dabs(grad_squared).gt.1.d-12)then
count_n += 1.d0
accu_relat += 2.0d0 * contrib/dabs(grad_squared+new)
endif
if(contrib.gt.1.d-10)then
print*,i,j,ipoint
print*,grad_squared,new,contrib
print*,2.0d0*contrib/dabs(grad_squared+new+1.d-12)
endif
accu += contrib
enddo
enddo
enddo
print*,'accu = ',accu/count_n
print*,'accu/rel = ',accu_relat/count_n
end
subroutine routine_fit
implicit none
integer :: i,nx
double precision :: dx,xmax,x,j_mu,j_mu_F_x_j,j_mu_fit_gauss
nx = 500
xmax = 5.d0
dx = xmax/dble(nx)
x = 0.d0
print*,'coucou',mu_erf
do i = 1, nx
write(33,'(100(F16.10,X))') x,j_mu(x),j_mu_F_x_j(x),j_mu_fit_gauss(x)
x += dx
enddo
end