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

Merge branch 'dev-stable-tc-scf' of github.com:QuantumPackage/qp2 into dev-stable-tc-scf

Browse Source
This commit is contained in:
Anthony Scemama 2023-03-28 13:18:22 +02:00
commit 405614c400
11 changed files with 277 additions and 169 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8

43
scripts/Hn.py Normal file
View File

@ -0,0 +1,43 @@
#!/usr/bin/env python
import sys
from math import *
arg = sys.argv
#f = open('data_dft','r')
n = int(sys.argv[1])
r = float(sys.argv[2])
f = open('H'+str(n)+'_'+str(r)+'.xyz','w')
string=str(n)+"\n"
f.write(string)
string="\n"
f.write(string)
for i in range(n):
x = r * cos(2.* i* pi/n)
y = r * sin(2.* i* pi/n)
z = 0.
string="H "+str(x)+" "+str(y)+" "+str(z)+"\n"
f.write(string)
#lines = f.readlines()
#cipsi_dft= []
#
#dissoc = []
#dissoc.append(float(-76.0179223470363))
#dissoc.append(float(-76.0592367866993))
#dissoc.append(float(-76.0678739715659))
#delta_e = []
#
#for line in lines:
# data = line.split()
# if(len(data)>0):
# dft=float(data[1])
# fci=float(data[2])
# e=fci+dft
# cipsi_dft.append(e)
#
#print(*cipsi_dft,sep=" & ")
#
#for i in 0,1,2:
# delta_e.append(1000.*(dissoc[i] - cipsi_dft[i]))
#
#print(*delta_e,sep=" & ")
#

2
scripts/get_fci_tc_conv.sh Executable file
View File

@ -0,0 +1,2 @@
file=$1
grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc

48
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f
View File

@ -38,7 +38,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, & !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, &
!$OMP ao_overlap_abs,sq_pi_3_2) !$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -46,7 +46,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
r(3) = final_grid_points(3,ipoint) r(3) = final_grid_points(3,ipoint)
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(ao_overlap_abs(j,i) .lt. 1.d-12) then if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
cycle cycle
endif endif
@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit) expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit)
if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle ! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j) int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
enddo enddo
@ -81,8 +81,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version ! if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
! expo_fit, i, j, int_fit_v, n_points_final_grid) ! expo_fit, i, j, int_fit_v, n_points_final_grid)
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -145,14 +144,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, big_array,& !$OMP List_comb_thr_b3_cent, big_array,&
!$OMP ao_abs_comb_b3_j1b,ao_overlap_abs) !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc)
! !
allocate(int_fit_v(n_points_final_grid)) allocate(int_fit_v(n_points_final_grid))
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(ao_overlap_abs(j,i) .lt. 1.d-12) then if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
cycle cycle
endif endif
@ -161,7 +160,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -243,7 +241,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, &
!$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b) !$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -260,11 +258,11 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i) int_j1b = ao_abs_comb_b3_j1b(1,j,i)
if(dabs(int_j1b).lt.1.d-10) cycle if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x_2(i_fit) expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit)
if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle ! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += coef_fit * int_fit tmp += coef_fit * int_fit
enddo enddo
@ -278,7 +276,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -288,8 +286,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef_fit = coef_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version ! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit tmp += coef * coef_fit * int_fit
enddo enddo
@ -350,7 +347,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2) !$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -369,7 +366,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -392,8 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv) sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
! if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version ! if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. 1d-10) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
@ -470,13 +466,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, & !$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test) !$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc) cycle
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint) r(2) = final_grid_points(2,ipoint)
@ -489,10 +485,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i) int_j1b = ao_abs_comb_b3_j1b(1,j,i)
if(dabs(int_j1b).lt.1.d-10) cycle ! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit) expo_fit = expo_gauss_j_mu_1_erf(i_fit)
if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle ! if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r)
tmp += coef_fit * int_fit tmp += coef_fit * int_fit
@ -507,7 +503,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -517,7 +513,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit) expo_fit = expo_gauss_j_mu_1_erf(i_fit)
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle ! if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)
alpha_1s = beta + expo_fit alpha_1s = beta + expo_fit
@ -527,9 +523,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3)) centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
if(expo_coef_1s .gt. 20.d0) cycle ! if(expo_coef_1s .gt. 20.d0) cycle
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
if(dabs(coef_tmp) .lt. 1d-08) cycle ! if(dabs(coef_tmp) .lt. 1d-08) cycle
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r) int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r)

46
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
View File

@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, &
!$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, & !$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
!do ipoint = 1, 10 !do ipoint = 1, 10
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
tmp = 0.d0 tmp = 0.d0
do i_1s = 1, List_comb_thr_b2_size(j,i) do i_1s = 1, List_comb_thr_b2_size(j,i)
@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -110,7 +110,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,& !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b, & !$OMP x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc)
! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss) ! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -120,7 +120,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
tmp_x = 0.d0 tmp_x = 0.d0
tmp_y = 0.d0 tmp_y = 0.d0
@ -130,19 +130,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
! if(ao_prod_center(1,j,i).ne.10000.d0)then
! ! approximate 1 - erf(mu r12) by a gaussian * 10
! !DIR$ FORCEINLINE
! call gaussian_product(expo_erfc_mu_gauss,r, &
! ao_prod_sigma(j,i),ao_prod_center(1,j,i), &
! factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-1.5d0)).lt.1.d-10)cycle
! endif
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
@ -216,7 +208,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -225,7 +217,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
tmp = 0.d0 tmp = 0.d0
@ -234,11 +226,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i) int_j1b = ao_abs_comb_b2_j1b(1,j,i)
if(dabs(int_j1b).lt.1.d-10) cycle ! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x(i_fit) expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit)
if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle ! if(ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc) cycle
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += coef_fit * int_fit tmp += coef_fit * int_fit
enddo enddo
@ -251,7 +243,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -259,9 +251,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
expo_fit = expo_gauss_j_mu_x(i_fit) expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit)
coeftot = coef * coef_fit coeftot = coef * coef_fit
if(dabs(coeftot).lt.1.d-15)cycle ! if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle ! if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit tmp += coef * coef_fit * int_fit
enddo enddo
@ -325,7 +317,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -334,7 +326,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
tmp = 0.d0 tmp = 0.d0
@ -343,7 +335,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i) int_j1b = ao_abs_comb_b2_j1b(1,j,i)
if(dabs(int_j1b).lt.1.d-10) cycle ! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
expo_fit = expo_good_j_mu_1gauss expo_fit = expo_good_j_mu_1gauss
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += int_fit tmp += int_fit
@ -356,7 +348,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -364,9 +356,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
expo_fit = expo_good_j_mu_1gauss expo_fit = expo_good_j_mu_1gauss
coef_fit = 1.d0 coef_fit = 1.d0
coeftot = coef * coef_fit coeftot = coef * coef_fit
if(dabs(coeftot).lt.1.d-15)cycle if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit tmp += coef * coef_fit * int_fit
! enddo ! enddo

42
src/ao_many_one_e_ints/listj1b_sorted.irp.f
View File

@ -3,15 +3,16 @@
&BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size] &BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size]
implicit none implicit none
integer :: i_1s,i,j,ipoint integer :: i_1s,i,j,ipoint
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-15
List_comb_thr_b2_size = 0 List_comb_thr_b2_size = 0
print*,'List_all_comb_b2_size = ',List_all_comb_b2_size
! pause
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
do i_1s = 1, List_all_comb_b2_size do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s) coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef).lt.1.d-15)cycle if(dabs(coef).lt.thrsh_cycle_tc)cycle
beta = List_all_comb_b2_expo (i_1s) beta = List_all_comb_b2_expo (i_1s)
beta = max(beta,1.d-12) beta = max(beta,1.d-12)
center(1:3) = List_all_comb_b2_cent(1:3,i_1s) center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
@ -24,7 +25,7 @@
dist += ( center(3) - r(3) )*( center(3) - r(3) ) dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo enddo
if(dabs(coef)*dabs(int_j1b).gt.thr)then if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
List_comb_thr_b2_size(j,i) += 1 List_comb_thr_b2_size(j,i) += 1
endif endif
enddo enddo
@ -40,6 +41,7 @@
list(i) = maxval(List_comb_thr_b2_size(:,i)) list(i) = maxval(List_comb_thr_b2_size(:,i))
enddo enddo
max_List_comb_thr_b2_size = maxval(list) max_List_comb_thr_b2_size = maxval(list)
print*,'max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size
END_PROVIDER END_PROVIDER
@ -49,16 +51,15 @@ END_PROVIDER
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)] &BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)]
implicit none implicit none
integer :: i_1s,i,j,ipoint,icount integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-15
ao_abs_comb_b2_j1b = 10000000.d0 ao_abs_comb_b2_j1b = 10000000.d0
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
icount = 0 icount = 0
do i_1s = 1, List_all_comb_b2_size do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s) coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef).lt.1.d-12)cycle if(dabs(coef).lt.thrsh_cycle_tc)cycle
beta = List_all_comb_b2_expo (i_1s) beta = List_all_comb_b2_expo (i_1s)
center(1:3) = List_all_comb_b2_cent(1:3,i_1s) center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
int_j1b = 0.d0 int_j1b = 0.d0
@ -70,7 +71,7 @@ END_PROVIDER
dist += ( center(3) - r(3) )*( center(3) - r(3) ) dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo enddo
if(dabs(coef)*dabs(int_j1b).gt.thr)then if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
icount += 1 icount += 1
List_comb_thr_b2_coef(icount,j,i) = coef List_comb_thr_b2_coef(icount,j,i) = coef
List_comb_thr_b2_expo(icount,j,i) = beta List_comb_thr_b2_expo(icount,j,i) = beta
@ -98,17 +99,17 @@ END_PROVIDER
&BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size] &BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size]
implicit none implicit none
integer :: i_1s,i,j,ipoint integer :: i_1s,i,j,ipoint
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-15
List_comb_thr_b3_size = 0 List_comb_thr_b3_size = 0
print*,'List_all_comb_b3_size = ',List_all_comb_b3_size
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
do i_1s = 1, List_all_comb_b3_size do i_1s = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s) coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s) beta = List_all_comb_b3_expo (i_1s)
center(1:3) = List_all_comb_b3_cent(1:3,i_1s) center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
if(dabs(coef).lt.thr)cycle if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_j1b = 0.d0 int_j1b = 0.d0
do ipoint = 1, n_points_extra_final_grid do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint) r(1:3) = final_grid_points_extra(1:3,ipoint)
@ -118,7 +119,7 @@ END_PROVIDER
dist += ( center(3) - r(3) )*( center(3) - r(3) ) dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo enddo
if(dabs(coef)*dabs(int_j1b).gt.thr)then if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
List_comb_thr_b3_size(j,i) += 1 List_comb_thr_b3_size(j,i) += 1
endif endif
enddo enddo
@ -144,9 +145,8 @@ END_PROVIDER
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)] &BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)]
implicit none implicit none
integer :: i_1s,i,j,ipoint,icount integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-15
ao_abs_comb_b3_j1b = 10000000.d0 ao_abs_comb_b3_j1b = 10000000.d0
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
@ -156,7 +156,7 @@ END_PROVIDER
beta = List_all_comb_b3_expo (i_1s) beta = List_all_comb_b3_expo (i_1s)
beta = max(beta,1.d-12) beta = max(beta,1.d-12)
center(1:3) = List_all_comb_b3_cent(1:3,i_1s) center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
if(dabs(coef).lt.thr)cycle if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_j1b = 0.d0 int_j1b = 0.d0
do ipoint = 1, n_points_extra_final_grid do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint) r(1:3) = final_grid_points_extra(1:3,ipoint)
@ -166,7 +166,7 @@ END_PROVIDER
dist += ( center(3) - r(3) )*( center(3) - r(3) ) dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo enddo
if(dabs(coef)*dabs(int_j1b).gt.thr)then if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
icount += 1 icount += 1
List_comb_thr_b3_coef(icount,j,i) = coef List_comb_thr_b3_coef(icount,j,i) = coef
List_comb_thr_b3_expo(icount,j,i) = beta List_comb_thr_b3_expo(icount,j,i) = beta
@ -177,15 +177,5 @@ END_PROVIDER
enddo enddo
enddo enddo
! do i = 1, ao_num
! do j = 1, i-1
! do icount = 1, List_comb_thr_b3_size(j,i)
! List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j)
! List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j)
! List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j)
! enddo
! enddo
! enddo
END_PROVIDER END_PROVIDER

19
src/non_h_ints_mu/total_tc_int.irp.f
View File

@ -68,7 +68,26 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)]
! --- ! ---
implicit none
integer :: i, j, k, l
double precision :: wall1, wall0
print *, ' providing ao_tc_int_chemist_no_cycle ...'
call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
! ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)] BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)]

6
src/tc_keywords/EZFIO.cfg
View File

@ -46,6 +46,12 @@ doc: Thresholds on the energy for iterative Davidson used in TC
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1.e-5 default: 1.e-5
[thrsh_cycle_tc]
type: Threshold
doc: Thresholds to cycle the integrals with the envelop
interface: ezfio,provider,ocaml
default: 1.e-15
[max_it_dav] [max_it_dav]
type: integer type: integer
doc: nb max of iteration in Davidson used in TC doc: nb max of iteration in Davidson used in TC

4
src/tc_scf/tc_scf.irp.f
View File

@ -11,8 +11,8 @@ program tc_scf
print *, ' starting ...' print *, ' starting ...'
my_grid_becke = .True. my_grid_becke = .True.
my_n_pt_r_grid = 30 my_n_pt_r_grid = 60
my_n_pt_a_grid = 50 my_n_pt_a_grid = 110
! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_r_grid = 10 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! 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 touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid

214
src/tc_scf/test_int.irp.f
View File

@ -21,25 +21,22 @@ program test_ints
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
!! OK !! OK
!call routine_int2_u_grad1u_j1b2 ! call routine_int2_u_grad1u_j1b2
!! OK ! OK
!call routine_v_ij_erf_rk_cst_mu_j1b ! call routine_v_ij_erf_rk_cst_mu_j1b
!! OK ! OK
! call routine_x_v_ij_erf_rk_cst_mu_j1b ! call routine_x_v_ij_erf_rk_cst_mu_j1b
!! OK ! OK
! call routine_v_ij_u_cst_mu_j1b ! call routine_int2_u2_j1b2
! OK
!! OK ! call routine_int2_u_grad1u_x_j1b2
!call routine_int2_u2_j1b2 ! OK
!! OK
!call routine_int2_u_grad1u_x_j1b2
!! OK
! call routine_int2_grad1u2_grad2u2_j1b2 ! call routine_int2_grad1u2_grad2u2_j1b2
! call routine_int2_u_grad1u_j1b2 ! call routine_int2_u_grad1u_j1b2
! call test_total_grad_lapl ! call test_total_grad_lapl
! call test_total_grad_square ! call test_total_grad_square
! call test_int2_grad1_u12_ao_test
! call routine_v_ij_u_cst_mu_j1b_test
! call test_ao_tc_int_chemist ! call test_ao_tc_int_chemist
! call test_grid_points_ao ! call test_grid_points_ao
! call test_tc_scf ! call test_tc_scf
@ -53,12 +50,12 @@ program test_ints
!call test_two_e_tc_non_hermit_integral() !call test_two_e_tc_non_hermit_integral()
call test_tc_grad_square_ao_test() ! call test_tc_grad_square_ao_test()
PROVIDE TC_HF_energy VARTC_HF_energy
print *, ' TC_HF_energy = ', TC_HF_energy
print *, ' VARTC_HF_energy = ', VARTC_HF_energy
!!PROVIDE TC_HF_energy VARTC_HF_energy
!!print *, ' TC_HF_energy = ', TC_HF_energy
!!print *, ' VARTC_HF_energy = ', VARTC_HF_energy
call test_old_ints
end end
! --- ! ---
@ -157,6 +154,9 @@ subroutine routine_int2_u_grad1u_j1b2
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_int2_u_grad1u_j1b2'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -169,20 +169,6 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b
integer :: i,j,ipoint,k,l integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
! print*,'ao_overlap_abs = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
! enddo
! print*,'center = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
! enddo
! print*,'sigma = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
! enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@ -215,6 +201,9 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_v_ij_erf_rk_cst_mu_j1b'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -228,20 +217,6 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b
integer :: i,j,ipoint,k,l,m integer :: i,j,ipoint,k,l,m
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
! print*,'ao_overlap_abs = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
! enddo
! print*,'center = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
! enddo
! print*,'sigma = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
! enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@ -276,6 +251,10 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_x_v_ij_erf_rk_cst_mu_j1b'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -290,20 +269,6 @@ subroutine routine_v_ij_u_cst_mu_j1b_test
integer :: i,j,ipoint,k,l integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
! print*,'ao_overlap_abs = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
! enddo
! print*,'center = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
! enddo
! print*,'sigma = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
! enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@ -336,6 +301,9 @@ subroutine routine_v_ij_u_cst_mu_j1b_test
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_v_ij_u_cst_mu_j1b_test'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -427,19 +395,6 @@ subroutine routine_int2_u2_j1b2
integer :: i,j,ipoint,k,l integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
! print*,'ao_overlap_abs = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
! enddo
! print*,'center = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
! enddo
! print*,'sigma = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
! enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
@ -473,6 +428,9 @@ subroutine routine_int2_u2_j1b2
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_int2_u2_j1b2'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -486,19 +444,6 @@ subroutine routine_int2_u_grad1u_x_j1b2
integer :: i,j,ipoint,k,l,m integer :: i,j,ipoint,k,l,m
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
! print*,'ao_overlap_abs = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
! enddo
! print*,'center = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
! enddo
! print*,'sigma = '
! do i = 1, ao_num
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
! enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
@ -534,6 +479,9 @@ subroutine routine_int2_u_grad1u_x_j1b2
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_int2_u_grad1u_x_j1b2'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -579,6 +527,9 @@ subroutine routine_v_ij_u_cst_mu_j1b
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'routine_v_ij_u_cst_mu_j1b'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -742,6 +693,9 @@ subroutine test_total_grad_lapl
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,' test_total_grad_lapl'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -767,6 +721,9 @@ subroutine test_total_grad_square
enddo enddo
enddo enddo
enddo enddo
print*,'******'
print*,'******'
print*,'test_total_grad_square'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4
@ -1057,3 +1014,86 @@ end
! --- ! ---
subroutine test_old_ints
implicit none
integer :: i,j,k,l
double precision :: old, new, contrib, get_ao_tc_sym_two_e_pot
double precision :: integral_sym , integral_nsym,accu
PROVIDE ao_tc_sym_two_e_pot_in_map
accu = 0.d0
do j = 1, ao_num
do l= 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
! integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
! old = integral_sym + integral_nsym
! old = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
new = ao_tc_int_chemist_test(k,i,l,j)
old = ao_tc_int_chemist_no_cycle(k,i,l,j)
contrib = dabs(old - new)
if(contrib.gt.1.d-6)then
print*,'problem !!'
print*,i,j,k,l
print*,old, new, contrib
endif
accu += contrib
enddo
enddo
enddo
enddo
print*,'******'
print*,'******'
print*,'in test_old_ints'
print*,'accu = ',accu/dble(ao_num**4)
end
subroutine test_int2_grad1_u12_ao_test
implicit none
integer :: i,j,ipoint,m,k,l
double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0
allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
array_ref = 0.d0
do m = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
do k = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do j = 1, ao_num
array(j,i,l,k) += int2_grad1_u12_ao_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
array_ref(j,i,l,k) += int2_grad1_u12_ao(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
enddo
enddo
enddo
enddo
enddo
enddo
accu_relat = 0.d0
accu_abs = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do j = 1, ao_num
contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
accu_abs += contrib
if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
accu_relat += contrib/dabs(array_ref(j,i,l,k))
endif
enddo
enddo
enddo
enddo
print*,'******'
print*,'******'
print*,'test_int2_grad1_u12_ao_test'
print*,'accu_abs = ',accu_abs/dble(ao_num)**4
print*,'accu_relat = ',accu_relat/dble(ao_num)**4
end

20
src/tools/sort_wf.irp.f Normal file
View File

@ -0,0 +1,20 @@
program sort_wf
implicit none
read_wf = .true.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.wf_sorted'
i_unit_output = getUnitAndOpen(output,'w')
do i = 1, N_det
write(i_unit_output, *)i,dabs(psi_coef_sorted(i,1))/dabs(psi_coef_sorted(1,1))
enddo
end