diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index 134cbcd2..3d7fbe50 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -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) ! 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 ! 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 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) BEGIN_DOC diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 088181cc..b7fe234f 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -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 coef_fit, expo_fit, int_fit, tmp) & !$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 List_all_comb_b3_coef, List_all_comb_b3_expo, & !$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 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) 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 coef_fit, expo_fit, int_fit, tmp) & !$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 List_all_comb_b3_coef, List_all_comb_b3_expo, & !$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 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) 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 tmp_x, tmp_y, tmp_z) & !$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 List_all_comb_b3_coef, List_all_comb_b3_expo, & !$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_y = 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) 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 alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & !$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 List_all_comb_b3_coef, List_all_comb_b3_expo, & !$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) 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) coef_fit = coef_gauss_j_mu_1_erf(i_fit) diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index 9486f153..6a662533 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -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 coef_fit, expo_fit, int_fit, tmp) & !$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 List_all_comb_b2_coef, List_all_comb_b2_expo, & !$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 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) coef_fit = coef_gauss_j_mu_x(i_fit) diff --git a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f index 5218f040..dddf98d4 100644 --- a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f @@ -575,17 +575,17 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A accu = 0.d0 ASSERT (n_pt_in > 1) - R1x(0) = (P_center(1) - A_center(1)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(1) - C_center(1))* p_new + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + 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 - R1xp(0) = (P_center(1) - B_center(1)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(1) - C_center(1))* p_new + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + 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 - R2x(0) = p_inv_2 - R2x(1) = 0.d0 - R2x(2) = -p_inv_2 * p_new + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2 * p_new !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 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 endif - R1x(0) = (P_center(2) - A_center(2)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(2) - C_center(2))* p_new + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + 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 - R1xp(0) = (P_center(2) - B_center(2)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(2) - C_center(2))* p_new + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + 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 a_y = power_A(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 endif - R1x(0) = (P_center(3) - A_center(3)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(3) - C_center(3)) * p_new + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + 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 - R1xp(0) = (P_center(3) - B_center(3)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(3) - C_center(3)) * p_new + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3)) * p_new !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 a_z = power_A(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 ! Returns the explicit polynomial in terms of the $t$ variable of the ! 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)$. END_DOC + implicit none integer, intent(in) :: n_pt_in integer,intent(out) :: n_pt_out diff --git a/src/ao_tc_eff_map/fit_j.irp.f b/src/ao_tc_eff_map/fit_j.irp.f index d53209d0..84ae28a7 100644 --- a/src/ao_tc_eff_map/fit_j.irp.f +++ b/src/ao_tc_eff_map/fit_j.irp.f @@ -1,3 +1,5 @@ +! --- + BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ] implicit none 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, coef_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, (ng_fit_jast)] BEGIN_DOC ! @@ -34,31 +36,71 @@ END_PROVIDER implicit none integer :: i 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 - call expo_fit_slater_gam(alpha, expos) - beta = expo_j_xmu(2) * mu_erf * mu_erf + coef_gauss_j_mu_x = (/ -0.47947881d0 /) + expo_gauss_j_mu_x = (/ 3.4987848d0 /) + + 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 - expo_gauss_j_mu_x(i) = expos(i) + beta - coef_gauss_j_mu_x(i) = tmp * coef_fit_slat_gauss(i) + endif + + 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 END_PROVIDER ! --- - BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x_2, (n_max_fit_slat)] -&BEGIN_PROVIDER [double precision, coef_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, (ng_fit_jast)] BEGIN_DOC ! ! 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 ! @@ -69,24 +111,64 @@ END_PROVIDER implicit none integer :: i 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 - !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 + if(ng_fit_jast .eq. 1) then - 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 - call expo_fit_slater_gam(alpha, expos) - beta = beta_opt * mu_erf * mu_erf + elseif(ng_fit_jast .eq. 2) then + + 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 - expo_gauss_j_mu_x_2(i) = expos(i) + beta + alpha = alpha_opt * mu_erf + 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) enddo diff --git a/src/ao_tc_eff_map/potential.irp.f b/src/ao_tc_eff_map/potential.irp.f index d8b3c442..8739b6cb 100644 --- a/src/ao_tc_eff_map/potential.irp.f +++ b/src/ao_tc_eff_map/potential.irp.f @@ -142,27 +142,72 @@ double precision function fit_1_erf_x(x) 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_DOC -! (1 - erf(mu*x))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * x^2) -! -! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2) -! -! and the slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians - END_DOC - integer :: i - double precision :: expos(n_max_fit_slat),alpha,beta - 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**2.d0 - 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 +! --- + + BEGIN_PROVIDER [double precision, expo_gauss_1_erf_x_2, (ng_fit_jast)] +&BEGIN_PROVIDER [double precision, coef_gauss_1_erf_x_2, (ng_fit_jast)] + + BEGIN_DOC + ! (1 - erf(mu*x))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * x^2) + ! + ! This is based on a fit of (1 - erf(mu*x)) by exp(-alpha * x) exp(-beta*mu^2x^2) + ! + ! and the slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians + END_DOC + + implicit none + integer :: i + double precision :: expos(ng_fit_jast), alpha, beta + + if(ng_fit_jast .eq. 1) then + + 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 +! --- + double precision function fit_1_erf_x_2(x) implicit none double precision, intent(in) :: x diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f index 0b8c78c2..926a20f1 100644 --- a/src/non_hermit_dav/biorthog.irp.f +++ b/src/non_hermit_dav/biorthog.irp.f @@ -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 n_good = 0 - thr = 1.d-5 + thr = 1.d-3 do i = 1, n print*, 'Re(i) + Im(i)', WR(i), WI(i) if(dabs(WI(i)) .lt. thr) then diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f index 4430109d..53c62ce8 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -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 @@ -1279,27 +1279,27 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap) integer :: i, j, num_linear_dependencies double precision :: threshold - double precision, allocatable :: S(:,:), tmp(:,:),Stmp(:,:) + double precision, allocatable :: S(:,:), tmp(:,:), Stmp(:,:) double precision, allocatable :: U(:,:), Vt(:,:), D(:) print *, ' apply SVD to orthogonalize vectors' ! --- - allocate(S(m,m),Stmp(n,m)) - ! 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) & , 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) & , 0.d0, S, size(S, 1) ) + deallocate(Stmp) -! print *, ' eigenvec overlap bef SVD: ' -! do i = 1, m -! write(*, '(1000(F16.10,X))') S(i,:) -! enddo + print *, ' eigenvec overlap bef SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo ! --- @@ -1348,23 +1348,23 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap) ! --- - allocate(S(m,m)) - - ! S = C.T x C - call dgemm( 'T', 'N', m, m, n, 1.d0 & - , C, size(C, 1), C, size(C, 1) & + ! S = C.T x overlap x C + allocate(S(m,m), Stmp(n,m)) + call dgemm( 'N', 'N', n, m, n, 1.d0 & + , overlap, size(overlap, 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) ) + deallocate(Stmp) -! print *, ' eigenvec overlap aft SVD: ' -! do i = 1, m -! write(*, '(1000(F16.10,X))') S(i,:) -! enddo - + print *, ' eigenvec overlap aft SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo deallocate(S) - ! --- - -end subroutine impose_orthog_svd +end subroutine impose_orthog_svd_overlap ! --- @@ -2780,8 +2780,6 @@ end subroutine check_weighted_biorthog_binormalize ! --- - - subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) implicit none @@ -2894,8 +2892,7 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) enddo deallocate(S) - ! --- - + return end subroutine impose_weighted_biorthog_svd ! --- diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 3cca5bb9..1afa26e9 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -118,7 +118,6 @@ doc: type of 1-body Jastrow interface: ezfio, provider, ocaml default: 0 - [thr_degen_tc] 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 @@ -130,3 +129,10 @@ type: logical doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors interface: ezfio,provider,ocaml default: False + +[ng_fit_jast] +type: integer +doc: nb of Gaussians used to fit Jastrow fcts +interface: ezfio,provider,ocaml +default: 6 + diff --git a/src/tc_scf/print_fit_param.irp.f b/src/tc_scf/print_fit_param.irp.f new file mode 100644 index 00000000..f8bcfa7f --- /dev/null +++ b/src/tc_scf/print_fit_param.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 1b03f338..42925e41 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -96,6 +96,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) ! n_degen = ilast - ifirst +1 n_degen = list_degen(i,0) + if(n_degen .eq. 1) cycle 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)) @@ -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(out) :: smat(n,n) integer :: i, j, k, l + double precision, allocatable :: S_tmp(:,:) - smat = 0.D0 - do i = 1, n - do j = 1, n - do k = 1, m - do l = 1, m - smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) - enddo - enddo - enddo - enddo + smat = 0.d0 + + !do i = 1, n + ! do j = 1, n + ! do k = 1, m + ! do l = 1, m + ! smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) + ! 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 ! --- -subroutine orthog_functions(m,n,coef,overlap) - implicit none - integer, intent(in) :: m,n - double precision, intent(in) :: overlap(m,m) - double precision, intent(inout) :: coef(m,n) - double precision, allocatable :: stmp(:,:) - integer :: j - allocate(stmp(n,n)) - call build_s_matrix(m,n,coef,coef,overlap,stmp) +subroutine orthog_functions(m, n, coef, overlap) + + implicit none + + integer, intent(in) :: m, n + double precision, intent(in) :: overlap(m,m) + double precision, intent(inout) :: coef(m,n) + double precision, allocatable :: stmp(:,:) + integer :: j + + allocate(stmp(n,n)) + call build_s_matrix(m, n, coef, coef, overlap, stmp) ! print*,'overlap before' ! do j = 1, n ! write(*,'(100(F16.10,X))')stmp(:,j) ! enddo - call impose_orthog_svd_overlap(m, n, coef,overlap) - call build_s_matrix(m,n,coef,coef,overlap,stmp) + call impose_orthog_svd_overlap(m, n, coef, overlap) + call build_s_matrix(m, n, coef, coef, overlap, stmp) do j = 1, n - coef(1,:m) *= 1.d0/dsqrt(stmp(j,j)) + coef(1,:m) *= 1.d0/dsqrt(stmp(j,j)) enddo - call build_s_matrix(m,n,coef,coef,overlap,stmp) + call build_s_matrix(m, n, coef, coef, overlap, stmp) !print*,'overlap after' !do j = 1, n diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 4a875b59..48cbbdc0 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -19,8 +19,8 @@ program tc_scf !call orthonormalize_mos call routine_scf() - call minimize_tc_orb_angles - call print_energy_and_mos + call minimize_tc_orb_angles() + call print_energy_and_mos() end diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f index 155eef9c..188bfa58 100644 --- a/src/utils/block_diag_degen.irp.f +++ b/src/utils/block_diag_degen.irp.f @@ -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 ! you enter with an array A(n) and spits out all the elements degenerated up to thr