LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-09 06:53:38 +01:00
Code Releases Activity

opt-g added

Browse Source
This commit is contained in:
AbdAmmar 2022-12-05 00:08:00 +01:00
parent a11de84bcd
commit d889e67146
13 changed files with 342 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ao_many_one_e_ints/ao_erf_gauss.irp.f
View File

@ -257,11 +257,11 @@ subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
! First term = (x-Ax)**(ax+1) ! First term = (x-Ax)**(ax+1)
integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in) integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in)
ints(m) += integral * coef ints(m) += integral * coef
! Second term = Ax * (x-Ax)**(ax) ! Second term = Ax * (x-Ax)**(ax)
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
ints(m) += A_center(m) * integral * coef ints(m) += A_center(m) * integral * coef
enddo enddo
@ -815,8 +815,6 @@ end function NAI_pol_x_mult_erf_ao_with1s_z
! --- ! ---
! ---
subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints) subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints)
BEGIN_DOC BEGIN_DOC

16
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -28,7 +28,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, & !$OMP final_grid_points, ng_fit_jast, &
!$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_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2) !$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
do j = i, ao_num do j = i, ao_num
tmp = 0.d0 tmp = 0.d0
do i_fit = 1, n_max_fit_slat 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 = coef_gauss_1_erf_x_2(i_fit) coef_fit = coef_gauss_1_erf_x_2(i_fit)
@ -120,7 +120,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, & !$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_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u2_j1b2) !$OMP List_all_comb_b3_cent, int2_u2_j1b2)
@ -134,7 +134,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
do j = i, ao_num do j = i, ao_num
tmp = 0.d0 tmp = 0.d0
do i_fit = 1, n_max_fit_slat 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)
@ -213,7 +213,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, & !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z) & !$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, & !$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_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2) !$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2)
@ -230,7 +230,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
tmp_x = 0.d0 tmp_x = 0.d0
tmp_y = 0.d0 tmp_y = 0.d0
tmp_z = 0.d0 tmp_z = 0.d0
do i_fit = 1, n_max_fit_slat 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)
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)
@ -330,7 +330,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, & !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, & !$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_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2) !$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2)
@ -343,7 +343,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
r(3) = final_grid_points(3,ipoint) r(3) = final_grid_points(3,ipoint)
tmp = 0.d0 tmp = 0.d0
do i_fit = 1, n_max_fit_slat 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)
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)

4
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -248,7 +248,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP final_grid_points, n_max_fit_slat, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b)
@ -263,7 +263,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
do j = i, ao_num do j = i, ao_num
tmp = 0.d0 tmp = 0.d0
do i_fit = 1, n_max_fit_slat 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)

47
src/ao_one_e_ints/pot_ao_erf_ints.irp.f
View File

@ -575,17 +575,17 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A
accu = 0.d0 accu = 0.d0
ASSERT (n_pt_in > 1) ASSERT (n_pt_in > 1)
R1x(0) = (P_center(1) - A_center(1)) R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0 R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))* p_new R1x(2) = -(P_center(1) - C_center(1))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(1) - B_center(1)) R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0 R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))* p_new R1xp(2) =-(P_center(1) - C_center(1))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R2x(0) = p_inv_2 R2x(0) = p_inv_2
R2x(1) = 0.d0 R2x(1) = 0.d0
R2x(2) = -p_inv_2 * p_new R2x(2) = -p_inv_2 * p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
do i = 0, n_pt_in do i = 0, n_pt_in
@ -609,13 +609,13 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A
return return
endif endif
R1x(0) = (P_center(2) - A_center(2)) R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0 R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))* p_new R1x(2) = -(P_center(2) - C_center(2))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(2) - B_center(2)) R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0 R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))* p_new R1xp(2) =-(P_center(2) - C_center(2))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
a_y = power_A(2) a_y = power_A(2)
b_y = power_B(2) b_y = power_B(2)
@ -628,13 +628,13 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A
return return
endif endif
R1x(0) = (P_center(3) - A_center(3)) R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0 R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3)) * p_new R1x(2) = -(P_center(3) - C_center(3)) * p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(3) - B_center(3)) R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0 R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3)) * p_new R1xp(2) =-(P_center(3) - C_center(3)) * p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
a_z = power_A(3) a_z = power_A(3)
b_z = power_B(3) b_z = power_B(3)
@ -663,16 +663,15 @@ end subroutine give_polynomial_mult_center_one_e_erf_opt
! --- ! ---
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,&
power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
BEGIN_DOC BEGIN_DOC
! Returns the explicit polynomial in terms of the $t$ variable of the ! Returns the explicit polynomial in terms of the $t$ variable of the
! following polynomial: ! following polynomial:
! !
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$. ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
END_DOC END_DOC
implicit none implicit none
integer, intent(in) :: n_pt_in integer, intent(in) :: n_pt_in
integer,intent(out) :: n_pt_out integer,intent(out) :: n_pt_out

134
src/ao_tc_eff_map/fit_j.irp.f
View File

@ -1,3 +1,5 @@
! ---
BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ] BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -14,8 +16,8 @@ END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x, (n_max_fit_slat)] BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x, (ng_fit_jast)]
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x, (n_max_fit_slat)] &BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x, (ng_fit_jast)]
BEGIN_DOC BEGIN_DOC
! !
@ -34,31 +36,71 @@ END_PROVIDER
implicit none implicit none
integer :: i integer :: i
double precision :: tmp double precision :: tmp
double precision :: expos(n_max_fit_slat), alpha, beta double precision :: expos(ng_fit_jast), alpha, beta
tmp = -0.5d0 / (mu_erf * sqrt(dacos(-1.d0))) if(ng_fit_jast .eq. 1) then
alpha = expo_j_xmu(1) * mu_erf coef_gauss_j_mu_x = (/ -0.47947881d0 /)
call expo_fit_slater_gam(alpha, expos) expo_gauss_j_mu_x = (/ 3.4987848d0 /)
beta = expo_j_xmu(2) * mu_erf * mu_erf
elseif(ng_fit_jast .eq. 2) then
coef_gauss_j_mu_x = (/ -0.18390742d0, -0.35512656d0 /)
expo_gauss_j_mu_x = (/ 31.9279947d0 , 2.11428789d0 /)
elseif(ng_fit_jast .eq. 3) then
coef_gauss_j_mu_x = (/ -0.07501725d0, -0.28499012d0, -0.1953932d0 /)
expo_gauss_j_mu_x = (/ 206.74058566d0, 1.72974157d0, 11.18735164d0 /)
elseif(ng_fit_jast .eq. 5) then
coef_gauss_j_mu_x = (/ -0.01832955d0 , -0.10188952d0 , -0.20710858d0 , -0.18975032d0 , -0.04641657d0 /)
expo_gauss_j_mu_x = (/ 4.33116687d+03, 2.61292842d+01, 1.43447161d+00, 4.92767426d+00, 2.10654699d+02 /)
elseif(ng_fit_jast .eq. 6) then
coef_gauss_j_mu_x = (/ -0.08783664d0 , -0.16088711d0 , -0.18464486d0 , -0.0368509d0 , -0.08130028d0 , -0.0126972d0 /)
expo_gauss_j_mu_x = (/ 4.09729729d+01, 7.11620618d+00, 2.03692338d+00, 4.10831731d+02, 1.12480198d+00, 1.00000000d+04 /)
elseif(ng_fit_jast .eq. 20) then
ASSERT(n_max_fit_slat == 20)
alpha = expo_j_xmu(1) * mu_erf
call expo_fit_slater_gam(alpha, expos)
beta = expo_j_xmu(2) * mu_erf * mu_erf
tmp = -1.0d0 / sqrt(dacos(-1.d0))
do i = 1, ng_fit_jast
expo_gauss_j_mu_x(i) = expos(i) + beta
coef_gauss_j_mu_x(i) = tmp * coef_fit_slat_gauss(i)
enddo
else
print *, ' not implemented yet'
stop
do i = 1, n_max_fit_slat endif
expo_gauss_j_mu_x(i) = expos(i) + beta
coef_gauss_j_mu_x(i) = tmp * coef_fit_slat_gauss(i) tmp = 0.5d0 / mu_erf
do i = 1, ng_fit_jast
coef_gauss_j_mu_x(i) = tmp * coef_gauss_j_mu_x(i)
enddo enddo
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x_2, (n_max_fit_slat)] BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x_2, (ng_fit_jast)]
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x_2, (n_max_fit_slat)] &BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x_2, (ng_fit_jast)]
BEGIN_DOC BEGIN_DOC
! !
! J(mu,r12)^2 = 0.25/mu^2 F(r12*mu)^2 ! J(mu,r12)^2 = 0.25/mu^2 F(r12*mu)^2
! !
! F(x)^2 = 1 /pi * exp(-2 * alpha * x) exp(-2 * beta * x^2) ! F(x)^2 = 1/pi * exp(-2 * alpha * x) exp(-2 * beta * x^2)
! !
! The slater function exp(-2 * alpha * x) is fitted with n_max_fit_slat gaussians ! The slater function exp(-2 * alpha * x) is fitted with n_max_fit_slat gaussians
! !
@ -69,24 +111,64 @@ END_PROVIDER
implicit none implicit none
integer :: i integer :: i
double precision :: tmp double precision :: tmp
double precision :: expos(n_max_fit_slat), alpha, beta double precision :: expos(ng_fit_jast), alpha, beta
double precision :: alpha_opt, beta_opt double precision :: alpha_opt, beta_opt
!alpha_opt = 2.d0 * expo_j_xmu(1) if(ng_fit_jast .eq. 1) then
!beta_opt = 2.d0 * expo_j_xmu(2)
! direct opt
alpha_opt = 3.52751759d0
beta_opt = 1.26214809d0
tmp = 0.25d0 / (mu_erf * mu_erf * dacos(-1.d0)) coef_gauss_j_mu_x_2 = (/ 0.26699573d0 /)
expo_gauss_j_mu_x_2 = (/ 11.71029824d0 /)
alpha = alpha_opt * mu_erf elseif(ng_fit_jast .eq. 2) then
call expo_fit_slater_gam(alpha, expos)
beta = beta_opt * mu_erf * mu_erf coef_gauss_j_mu_x_2 = (/ 0.11627934d0 , 0.18708824d0 /)
expo_gauss_j_mu_x_2 = (/ 102.41386863d0, 6.36239771d0 /)
elseif(ng_fit_jast .eq. 3) then
coef_gauss_j_mu_x_2 = (/ 0.04947216d0 , 0.14116238d0, 0.12276501d0 /)
expo_gauss_j_mu_x_2 = (/ 635.29701766d0, 4.87696954d0, 33.36745891d0 /)
elseif(ng_fit_jast .eq. 5) then
coef_gauss_j_mu_x_2 = (/ 0.01461527d0 , 0.03257147d0 , 0.08831354d0 , 0.11411794d0 , 0.06858783d0 /)
expo_gauss_j_mu_x_2 = (/ 8.76554470d+03, 4.90224577d+02, 3.68267125d+00, 1.29663940d+01, 6.58240931d+01 /)
elseif(ng_fit_jast .eq. 6) then
coef_gauss_j_mu_x_2 = (/ 0.01347632d0 , 0.03929124d0 , 0.06289468d0 , 0.10702493d0 , 0.06999865d0 , 0.02558191d0 /)
expo_gauss_j_mu_x_2 = (/ 1.00000000d+04, 1.20900717d+02, 3.20346191d+00, 8.92157196d+00, 3.28119120d+01, 6.49045808d+02 /)
elseif(ng_fit_jast .eq. 20) then
ASSERT(n_max_fit_slat == 20)
!alpha_opt = 2.d0 * expo_j_xmu(1)
!beta_opt = 2.d0 * expo_j_xmu(2)
! direct opt
alpha_opt = 3.52751759d0
beta_opt = 1.26214809d0
do i = 1, n_max_fit_slat alpha = alpha_opt * mu_erf
expo_gauss_j_mu_x_2(i) = expos(i) + beta call expo_fit_slater_gam(alpha, expos)
beta = beta_opt * mu_erf * mu_erf
tmp = 1.d0 / dacos(-1.d0)
do i = 1, ng_fit_jast
expo_gauss_j_mu_x_2(i) = expos(i) + beta
coef_gauss_j_mu_x_2(i) = tmp * coef_fit_slat_gauss(i)
enddo
else
print *, ' not implemented yet'
stop
endif
tmp = 0.25d0 / (mu_erf * mu_erf)
do i = 1, ng_fit_jast
coef_gauss_j_mu_x_2(i) = tmp * coef_fit_slat_gauss(i) coef_gauss_j_mu_x_2(i) = tmp * coef_fit_slat_gauss(i)
enddo enddo

83
src/ao_tc_eff_map/potential.irp.f
View File

@ -142,27 +142,72 @@ double precision function fit_1_erf_x(x)
end end
BEGIN_PROVIDER [double precision, expo_gauss_1_erf_x_2, (n_max_fit_slat)] ! ---
&BEGIN_PROVIDER [double precision, coef_gauss_1_erf_x_2, (n_max_fit_slat)]
implicit none BEGIN_PROVIDER [double precision, expo_gauss_1_erf_x_2, (ng_fit_jast)]
BEGIN_DOC &BEGIN_PROVIDER [double precision, coef_gauss_1_erf_x_2, (ng_fit_jast)]
! (1 - erf(mu*x))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * x^2)
! BEGIN_DOC
! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2) ! (1 - erf(mu*x))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * x^2)
! !
! and the slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians ! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2)
END_DOC !
integer :: i ! and the slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians
double precision :: expos(n_max_fit_slat),alpha,beta END_DOC
alpha = 2.d0 * expos_slat_gauss_1_erf_x(1) * mu_erf
call expo_fit_slater_gam(alpha,expos) implicit none
beta = 2.d0 * expos_slat_gauss_1_erf_x(2) * mu_erf**2.d0 integer :: i
do i = 1, n_max_fit_slat double precision :: expos(ng_fit_jast), alpha, beta
expo_gauss_1_erf_x_2(i) = expos(i) + beta
coef_gauss_1_erf_x_2(i) = coef_fit_slat_gauss(i) if(ng_fit_jast .eq. 1) then
enddo
coef_gauss_j_mu_x = (/ 0.85345277d0 /)
expo_gauss_j_mu_x = (/ 6.23519457d0 /)
elseif(ng_fit_jast .eq. 2) then
coef_gauss_j_mu_x = (/ 0.31030624d0 , 0.64364964d0 /)
expo_gauss_j_mu_x = (/ 55.39184787d0, 3.92151407d0 /)
elseif(ng_fit_jast .eq. 3) then
coef_gauss_j_mu_x = (/ 0.33206082d0 , 0.52347449d0, 0.12605012d0 /)
expo_gauss_j_mu_x = (/ 19.90272209d0, 3.2671671d0 , 336.47320445d0 /)
elseif(ng_fit_jast .eq. 5) then
coef_gauss_j_mu_x = (/ 0.02956716d0, 0.17025555d0, 0.32774114d0, 0.39034764d0, 0.07822781d0 /)
expo_gauss_j_mu_x = (/ 6467.28126d0, 46.9071990d0, 9.09617721d0, 2.76883328d0, 360.367093d0 /)
elseif(ng_fit_jast .eq. 6) then
coef_gauss_j_mu_x = (/ 0.18331042d0 , 0.10971118d0 , 0.29949169d0 , 0.34853132d0 , 0.0394275d0 , 0.01874444d0 /)
expo_gauss_j_mu_x = (/ 2.54293498d+01, 1.40317872d+02, 7.14630801d+00, 2.65517675d+00, 1.45142619d+03, 1.00000000d+04 /)
elseif(ng_fit_jast .eq. 20) then
ASSERT(n_max_fit_slat == 20)
alpha = 2.d0 * expos_slat_gauss_1_erf_x(1) * mu_erf
call expo_fit_slater_gam(alpha, expos)
beta = 2.d0 * expos_slat_gauss_1_erf_x(2) * mu_erf * mu_erf
do i = 1, n_max_fit_slat
expo_gauss_1_erf_x_2(i) = expos(i) + beta
coef_gauss_1_erf_x_2(i) = coef_fit_slat_gauss(i)
enddo
else
print *, ' not implemented yet'
stop
endif
END_PROVIDER END_PROVIDER
! ---
double precision function fit_1_erf_x_2(x) double precision function fit_1_erf_x_2(x)
implicit none implicit none
double precision, intent(in) :: x double precision, intent(in) :: x

2
src/non_hermit_dav/biorthog.irp.f
View File

@ -328,7 +328,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
! track & sort the real eigenvalues ! track & sort the real eigenvalues
n_good = 0 n_good = 0
thr = 1.d-5 thr = 1.d-3
do i = 1, n do i = 1, n
print*, 'Re(i) + Im(i)', WR(i), WI(i) print*, 'Re(i) + Im(i)', WR(i), WI(i)
if(dabs(WI(i)) .lt. thr) then if(dabs(WI(i)) .lt. thr) then

51
src/non_hermit_dav/lapack_diag_non_hermit.irp.f
View File

@ -1269,7 +1269,7 @@ end subroutine impose_orthog_svd
! --- ! ---
subroutine impose_orthog_svd_overlap(n, m, C,overlap) subroutine impose_orthog_svd_overlap(n, m, C, overlap)
implicit none implicit none
@ -1279,27 +1279,27 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
integer :: i, j, num_linear_dependencies integer :: i, j, num_linear_dependencies
double precision :: threshold double precision :: threshold
double precision, allocatable :: S(:,:), tmp(:,:),Stmp(:,:) double precision, allocatable :: S(:,:), tmp(:,:), Stmp(:,:)
double precision, allocatable :: U(:,:), Vt(:,:), D(:) double precision, allocatable :: U(:,:), Vt(:,:), D(:)
print *, ' apply SVD to orthogonalize vectors' print *, ' apply SVD to orthogonalize vectors'
! --- ! ---
allocate(S(m,m),Stmp(n,m))
! S = C.T x overlap x C ! S = C.T x overlap x C
call dgemm( 'N', 'N', n, m, n, 1.d0 & allocate(S(m,m), Stmp(n,m))
call dgemm( 'N', 'N', n, m, n, 1.d0 &
, overlap, size(overlap, 1), C, size(C, 1) & , overlap, size(overlap, 1), C, size(C, 1) &
, 0.d0, Stmp, size(Stmp, 1) ) , 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'T', 'N', m, m, n, 1.d0 & call dgemm( 'T', 'N', m, m, n, 1.d0 &
, C, size(C, 1), Stmp, size(Stmp, 1) & , C, size(C, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) ) , 0.d0, S, size(S, 1) )
deallocate(Stmp)
! print *, ' eigenvec overlap bef SVD: ' print *, ' eigenvec overlap bef SVD: '
! do i = 1, m do i = 1, m
! write(*, '(1000(F16.10,X))') S(i,:) write(*, '(1000(F16.10,X))') S(i,:)
! enddo enddo
! --- ! ---
@ -1348,23 +1348,23 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
! --- ! ---
allocate(S(m,m)) ! S = C.T x overlap x C
allocate(S(m,m), Stmp(n,m))
! S = C.T x C call dgemm( 'N', 'N', n, m, n, 1.d0 &
call dgemm( 'T', 'N', m, m, n, 1.d0 & , overlap, size(overlap, 1), C, size(C, 1) &
, C, size(C, 1), C, size(C, 1) & , 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, C, size(C, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) ) , 0.d0, S, size(S, 1) )
deallocate(Stmp)
! print *, ' eigenvec overlap aft SVD: ' print *, ' eigenvec overlap aft SVD: '
! do i = 1, m do i = 1, m
! write(*, '(1000(F16.10,X))') S(i,:) write(*, '(1000(F16.10,X))') S(i,:)
! enddo enddo
deallocate(S) deallocate(S)
! --- end subroutine impose_orthog_svd_overlap
end subroutine impose_orthog_svd
! --- ! ---
@ -2780,8 +2780,6 @@ end subroutine check_weighted_biorthog_binormalize
! --- ! ---
subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
implicit none implicit none
@ -2894,8 +2892,7 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
enddo enddo
deallocate(S) deallocate(S)
! --- return
end subroutine impose_weighted_biorthog_svd end subroutine impose_weighted_biorthog_svd
! --- ! ---

8
src/tc_keywords/EZFIO.cfg
View File

@ -118,7 +118,6 @@ doc: type of 1-body Jastrow
interface: ezfio, provider, ocaml interface: ezfio, provider, ocaml
default: 0 default: 0
[thr_degen_tc] [thr_degen_tc]
type: Threshold type: Threshold
doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue
@ -130,3 +129,10 @@ type: logical
doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[ng_fit_jast]
type: integer
doc: nb of Gaussians used to fit Jastrow fcts
interface: ezfio,provider,ocaml
default: 6

60
src/tc_scf/print_fit_param.irp.f Normal file
View File

@ -0,0 +1,60 @@
program print_fit_param
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
! 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 create_guess
!call orthonormalize_mos
call main()
end
! ---
subroutine main()
implicit none
integer :: i
mu_erf = 1.d0
touch mu_erf
print *, ' fit for (1 - erf(x))^2'
do i = 1, n_max_fit_slat
print*, expo_gauss_1_erf_x_2(i), coef_gauss_1_erf_x_2(i)
enddo
print *, ''
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]'
do i = 1, n_max_fit_slat
print *, expo_gauss_j_mu_x(i), 2.d0 * coef_gauss_j_mu_x(i)
enddo
print *, ''
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]^2'
do i = 1, n_max_fit_slat
print *, expo_gauss_j_mu_x_2(i), 4.d0 * coef_gauss_j_mu_x_2(i)
enddo
print *, ''
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)] x [1 - erf(mu * r12)]'
do i = 1, n_max_fit_slat
print *, expo_gauss_j_mu_1_erf(i), 4.d0 * coef_gauss_j_mu_1_erf(i)
enddo
return
end subroutine main
! ---

63
src/tc_scf/routines_rotates.irp.f
View File

@ -96,6 +96,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
! n_degen = ilast - ifirst +1 ! n_degen = ilast - ifirst +1
n_degen = list_degen(i,0) n_degen = list_degen(i,0)
if(n_degen .eq. 1) cycle
allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen)) allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen)) allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
@ -221,41 +222,57 @@ subroutine build_s_matrix(m, n, C1, C2, overlap, smat)
double precision, intent(in) :: C1(m,n), C2(m,n), overlap(m,m) double precision, intent(in) :: C1(m,n), C2(m,n), overlap(m,m)
double precision, intent(out) :: smat(n,n) double precision, intent(out) :: smat(n,n)
integer :: i, j, k, l integer :: i, j, k, l
double precision, allocatable :: S_tmp(:,:)
smat = 0.D0 smat = 0.d0
do i = 1, n
do j = 1, n !do i = 1, n
do k = 1, m ! do j = 1, n
do l = 1, m ! do k = 1, m
smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) ! do l = 1, m
enddo ! smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j)
enddo ! enddo
enddo ! enddo
enddo ! enddo
!enddo
! C1.T x overlap
allocate(S_tmp(n,m))
call dgemm( 'T', 'N', n, m, m, 1.d0 &
, C1, size(C1, 1), overlap, size(overlap, 1) &
, 0.d0, S_tmp, size(S_tmp, 1) )
! C1.T x overlap x C2
call dgemm( 'N', 'N', n, n, m, 1.d0 &
, S_tmp, size(S_tmp, 1), C2(1,1), size(C2, 1) &
, 0.d0, smat, size(smat, 1) )
deallocate(S_tmp)
end end
! --- ! ---
subroutine orthog_functions(m,n,coef,overlap) subroutine orthog_functions(m, n, coef, overlap)
implicit none
integer, intent(in) :: m,n implicit none
double precision, intent(in) :: overlap(m,m)
double precision, intent(inout) :: coef(m,n) integer, intent(in) :: m, n
double precision, allocatable :: stmp(:,:) double precision, intent(in) :: overlap(m,m)
integer :: j double precision, intent(inout) :: coef(m,n)
allocate(stmp(n,n)) double precision, allocatable :: stmp(:,:)
call build_s_matrix(m,n,coef,coef,overlap,stmp) integer :: j
allocate(stmp(n,n))
call build_s_matrix(m, n, coef, coef, overlap, stmp)
! print*,'overlap before' ! print*,'overlap before'
! do j = 1, n ! do j = 1, n
! write(*,'(100(F16.10,X))')stmp(:,j) ! write(*,'(100(F16.10,X))')stmp(:,j)
! enddo ! enddo
call impose_orthog_svd_overlap(m, n, coef,overlap) call impose_orthog_svd_overlap(m, n, coef, overlap)
call build_s_matrix(m,n,coef,coef,overlap,stmp) call build_s_matrix(m, n, coef, coef, overlap, stmp)
do j = 1, n do j = 1, n
coef(1,:m) *= 1.d0/dsqrt(stmp(j,j)) coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
enddo enddo
call build_s_matrix(m,n,coef,coef,overlap,stmp) call build_s_matrix(m, n, coef, coef, overlap, stmp)
!print*,'overlap after' !print*,'overlap after'
!do j = 1, n !do j = 1, n

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

@ -19,8 +19,8 @@ program tc_scf
!call orthonormalize_mos !call orthonormalize_mos
call routine_scf() call routine_scf()
call minimize_tc_orb_angles call minimize_tc_orb_angles()
call print_energy_and_mos call print_energy_and_mos()
end end

2
src/utils/block_diag_degen.irp.f
View File

@ -155,7 +155,7 @@ end
! --- ! ---
subroutine give_degen_full_list(A,n,thr,list_degen,n_degen_list) subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list)
BEGIN_DOC BEGIN_DOC
! you enter with an array A(n) and spits out all the elements degenerated up to thr ! you enter with an array A(n) and spits out all the elements degenerated up to thr