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 39be352f..d9c35a8c 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 @@ -1,4 +1,6 @@ +! --- + subroutine phi_j_erf_mu_r_xyz_phi(i,j,mu_in, C_center, xyz_ints) implicit none BEGIN_DOC @@ -49,45 +51,58 @@ subroutine phi_j_erf_mu_r_xyz_phi(i,j,mu_in, C_center, xyz_ints) enddo end +! --- -double precision function phi_j_erf_mu_r_phi(i,j,mu_in, C_center) - implicit none - BEGIN_DOC -! phi_j_erf_mu_r_phi = int dr phi_j(r) [erf(mu |r - C|)/|r-C|] phi_i(r) - END_DOC - integer, intent(in) :: i,j - double precision, intent(in) :: mu_in, C_center(3) - integer :: num_A,power_A(3), num_b, power_B(3) - double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf - integer :: n_pt_in,l,m - phi_j_erf_mu_r_phi = 0.d0 - if(ao_overlap_abs(j,i).lt.1.d-12)then - return - endif - n_pt_in = n_pt_max_integrals - ! j - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - ! i - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) +double precision function phi_j_erf_mu_r_phi(i, j, mu_in, C_center) - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - contrib = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in) - phi_j_erf_mu_r_phi += contrib * ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) + BEGIN_DOC + ! phi_j_erf_mu_r_phi = int dr phi_j(r) [erf(mu |r - C|)/|r-C|] phi_i(r) + END_DOC + + implicit none + integer, intent(in) :: i,j + double precision, intent(in) :: mu_in, C_center(3) + + integer :: num_A, power_A(3), num_b, power_B(3) + integer :: n_pt_in, l, m + double precision :: alpha, beta, A_center(3), B_center(3), contrib + + double precision :: NAI_pol_mult_erf + + phi_j_erf_mu_r_phi = 0.d0 + + if(ao_overlap_abs(j,i).lt.1.d-12) then + return + endif + + n_pt_in = n_pt_max_integrals + + ! j + num_A = ao_nucl(j) + power_A(1:3) = ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + ! i + num_B = ao_nucl(i) + power_B(1:3) = ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l = 1, ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + do m = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + + contrib = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + phi_j_erf_mu_r_phi += contrib * ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i) + enddo enddo - enddo -end +end function phi_j_erf_mu_r_phi +! --- -subroutine erfc_mu_gauss_xyz_ij_ao(i,j,mu, C_center, delta,gauss_ints) +subroutine erfc_mu_gauss_xyz_ij_ao(i, j, mu, C_center, delta, gauss_ints) implicit none BEGIN_DOC ! gauss_ints(m) = \int dr exp(-delta (r - C)^2 ) x/y/z * ( 1 - erf(mu |r-r'|))/ |r-r'| * AO_i(r') * AO_j(r') @@ -132,95 +147,204 @@ subroutine erfc_mu_gauss_xyz_ij_ao(i,j,mu, C_center, delta,gauss_ints) enddo end -subroutine erf_mu_gauss_ij_ao(i,j,mu, C_center, delta,gauss_ints) - implicit none - BEGIN_DOC - ! gauss_ints(m) = \int dr exp(-delta (r - C)^2 ) * erf(mu |r-r'|)/ |r-r'| * AO_i(r') * AO_j(r') - ! - END_DOC - integer, intent(in) :: i,j - double precision, intent(in) :: mu, C_center(3),delta - double precision, intent(out):: gauss_ints +! --- - integer :: num_A,power_A(3), num_b, power_B(3) - double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf - double precision :: integral , erf_mu_gauss - integer :: n_pt_in,l,m,mm - gauss_ints = 0.d0 - if(ao_overlap_abs(j,i).lt.1.d-12)then - return - endif - n_pt_in = n_pt_max_integrals - ! j - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - ! i - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) +subroutine erf_mu_gauss_ij_ao(i, j, mu, C_center, delta, gauss_ints) - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - if(dabs(ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i)).lt.1.d-12)cycle - integral = erf_mu_gauss(C_center,delta,mu,A_center,B_center,power_A,power_B,alpha,beta,n_pt_in) - gauss_ints += integral * ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) - enddo - enddo -end - - -subroutine NAI_pol_x_mult_erf_ao(i_ao,j_ao,mu_in,C_center,ints) - implicit none BEGIN_DOC + ! + ! gauss_ints = \int dr exp(-delta (r - C)^2) * erf(mu |r-C|) / |r-C| * AO_i(r) * AO_j(r) + ! + END_DOC + + implicit none + integer, intent(in) :: i, j + double precision, intent(in) :: mu, C_center(3), delta + double precision, intent(out) :: gauss_ints + + integer :: n_pt_in, l, m + integer :: num_A, power_A(3), num_b, power_B(3) + double precision :: alpha, beta, A_center(3), B_center(3), coef + double precision :: integral + + double precision :: erf_mu_gauss + + gauss_ints = 0.d0 + + if(ao_overlap_abs(j,i).lt.1.d-12) then + return + endif + + n_pt_in = n_pt_max_integrals + + ! j + num_A = ao_nucl(j) + power_A(1:3) = ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + ! i + num_B = ao_nucl(i) + power_B(1:3) = ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l = 1, ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + do m = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + coef = ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i) + + if(dabs(coef) .lt. 1.d-12) cycle + + integral = erf_mu_gauss(C_center, delta, mu, A_center, B_center, power_A, power_B, alpha, beta, n_pt_in) + + gauss_ints += integral * coef + enddo + enddo + +end subroutine erf_mu_gauss_ij_ao + +! --- + +subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + + BEGIN_DOC + ! ! Computes the following integral : + ! ! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! ! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! ! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! END_DOC - include 'utils/constants.include.F' - integer, intent(in) :: i_ao,j_ao - double precision, intent(in) :: mu_in, C_center(3) - double precision, intent(out):: ints(3) - double precision :: A_center(3), B_center(3),integral, alpha,beta - double precision :: NAI_pol_mult_erf - integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in, power_xA(3),m - ints = 0.d0 - if(ao_overlap_abs(j_ao,i_ao).lt.1.d-12)then - return - endif - num_A = ao_nucl(i_ao) - power_A(1:3)= ao_power(i_ao,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - num_B = ao_nucl(j_ao) - power_B(1:3)= ao_power(j_ao,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - n_pt_in = n_pt_max_integrals + include 'utils/constants.include.F' - do i = 1, ao_prim_num(i_ao) - alpha = ao_expo_ordered_transp(i,i_ao) - do m = 1, 3 - power_xA = power_A - ! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax - power_xA(m) += 1 - do j = 1, ao_prim_num(j_ao) - beta = ao_expo_ordered_transp(j,j_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) - ints(m) += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao) - ! 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) - ints(m) += A_center(m) * integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao) + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + double precision, intent(out) :: ints(3) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, power_xA(3), m + double precision :: A_center(3), B_center(3), integral, alpha, beta, coef + + double precision :: NAI_pol_mult_erf + + ints = 0.d0 + if(ao_overlap_abs(j_ao,i_ao).lt.1.d-12) then + return + endif + + num_A = ao_nucl(i_ao) + power_A(1:3) = ao_power(i_ao,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(j_ao) + power_B(1:3) = ao_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alpha = ao_expo_ordered_transp(i,i_ao) + + do m = 1, 3 + + power_xA = power_A + ! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax + power_xA(m) += 1 + + do j = 1, ao_prim_num(j_ao) + beta = ao_expo_ordered_transp(j,j_ao) + 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) + 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) + ints(m) += A_center(m) * integral * coef + + enddo enddo enddo - enddo -end + +end subroutine NAI_pol_x_mult_erf_ao + +! --- + +subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3) + double precision, intent(out) :: ints(3) + + integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3), m + double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef + + double precision, external :: NAI_pol_mult_erf_with1s + + ints = 0.d0 + if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then + return + endif + + power_Ai(1:3) = ao_power(i_ao,1:3) + power_Aj(1:3) = ao_power(j_ao,1:3) + + Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3) + Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alphai = ao_expo_ordered_transp(i,i_ao) + + do m = 1, 3 + + power_xA = power_Ai + ! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax + power_xA(m) += 1 + + do j = 1, ao_prim_num(j_ao) + alphaj = ao_expo_ordered_transp(j,j_ao) + 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_with1s( Ai_center, Aj_center, power_xA, power_Aj, alphai, alphaj & + , beta, b_center, c_center, n_pt_in, mu_in ) + ints(m) += integral * coef + + ! Second term = Ax * (x-Ax)**(ax) + integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj & + , beta, b_center, c_center, n_pt_in, mu_in ) + ints(m) += Ai_center(m) * integral * coef + + enddo + enddo + enddo + +end subroutine NAI_pol_x_mult_erf_ao_with1s + +! --- subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints) implicit none @@ -249,7 +373,6 @@ subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints) B_center(1:3) = nucl_coord(num_B,1:3) n_pt_in = n_pt_max_integrals - do i = 1, ao_prim_num(i_ao) alpha = ao_expo_ordered_transp(i,i_ao) power_xA = power_A @@ -267,3 +390,5 @@ subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints) enddo end +! --- + diff --git a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f index 681d1e6f..fadec343 100644 --- a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f @@ -102,36 +102,118 @@ subroutine overlap_gauss_r12_all_ao(D_center,delta,aos_ints) enddo end -double precision function overlap_gauss_r12_ao(D_center,delta,i,j) - implicit none - BEGIN_DOC -! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2} - END_DOC - integer, intent(in) :: i,j - double precision, intent(in) :: D_center(3), delta +! --- - integer :: num_a,num_b,power_A(3), power_B(3),l,k - double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,analytical_j - overlap_gauss_r12_ao = 0.d0 - if(ao_overlap_abs(j,i).lt.1.d-12)then - return - endif - ! TODO :: PUT CYCLES IN LOOPS - num_A = ao_nucl(i) - power_A(1:3)= ao_power(i,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - num_B = ao_nucl(j) - power_B(1:3)= ao_power(j,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - do l=1,ao_prim_num(i) - alpha = ao_expo_ordered_transp(l,i) - do k=1,ao_prim_num(j) - beta = ao_expo_ordered_transp(k,j) - analytical_j = overlap_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta) - overlap_gauss_r12_ao += analytical_j * ao_coef_normalized_ordered_transp(l,i) & - * ao_coef_normalized_ordered_transp(k,j) - enddo - enddo -end +! TODO :: PUT CYCLES IN LOOPS +double precision function overlap_gauss_r12_ao(D_center, delta, i, j) + BEGIN_DOC + ! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2} + END_DOC + + implicit none + integer, intent(in) :: i, j + double precision, intent(in) :: D_center(3), delta + + integer :: power_A(3), power_B(3), l, k + double precision :: A_center(3), B_center(3), alpha, beta, coef, analytical_j + + double precision, external :: overlap_gauss_r12 + + overlap_gauss_r12_ao = 0.d0 + + if(ao_overlap_abs(j,i).lt.1.d-12) then + return + endif + + power_A(1:3) = ao_power(i,1:3) + power_B(1:3) = ao_power(j,1:3) + + A_center(1:3) = nucl_coord(ao_nucl(i),1:3) + B_center(1:3) = nucl_coord(ao_nucl(j),1:3) + + do l = 1, ao_prim_num(i) + alpha = ao_expo_ordered_transp(l,i) + do k = 1, ao_prim_num(j) + beta = ao_expo_ordered_transp(k,j) + coef = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j) + + if(dabs(coef) .lt. 1d-12) cycle + + analytical_j = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta) + + overlap_gauss_r12_ao += coef * analytical_j + enddo + enddo + +end function overlap_gauss_r12_ao + +! --- + +double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j) + + BEGIN_DOC + ! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} + END_DOC + + implicit none + integer, intent(in) :: i, j + double precision, intent(in) :: B_center(3), beta, D_center(3), delta + + integer :: power_A1(3), power_A2(3), l, k + double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, analytical_j + double precision :: G_center(3), gama, fact_g, gama_inv + + double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao + + ASSERT(beta .gt. 0.d0) + if(beta .lt. 1d-10) then + overlap_gauss_r12_ao_with1s = overlap_gauss_r12_ao(D_center, delta, i, j) + return + endif + + overlap_gauss_r12_ao_with1s = 0.d0 + + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + return + endif + + ! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2} + + gama = beta + delta + gama_inv = 1.d0 / gama + G_center(1) = (beta * B_center(1) + delta * D_center(1)) * gama_inv + G_center(2) = (beta * B_center(2) + delta * D_center(2)) * gama_inv + G_center(3) = (beta * B_center(3) + delta * D_center(3)) * gama_inv + fact_g = beta * delta * gama_inv * ( (B_center(1) - D_center(1)) * (B_center(1) - D_center(1)) & + + (B_center(2) - D_center(2)) * (B_center(2) - D_center(2)) & + + (B_center(3) - D_center(3)) * (B_center(3) - D_center(3)) ) + fact_g = dexp(-fact_g) + if(fact_g .lt. 1.d-12) return + + ! --- + + power_A1(1:3) = ao_power(i,1:3) + power_A2(1:3) = ao_power(j,1:3) + + A1_center(1:3) = nucl_coord(ao_nucl(i),1:3) + A2_center(1:3) = nucl_coord(ao_nucl(j),1:3) + + do l = 1, ao_prim_num(i) + alpha1 = ao_expo_ordered_transp(l,i) + do k = 1, ao_prim_num(j) + alpha2 = ao_expo_ordered_transp(k,j) + coef12 = fact_g * ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j) + + if(dabs(coef12) .lt. 1d-12) cycle + + analytical_j = overlap_gauss_r12(G_center, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2) + + overlap_gauss_r12_ao_with1s += coef12 * analytical_j + enddo + enddo + +end function overlap_gauss_r12_ao_with1s + +! --- diff --git a/src/ao_many_one_e_ints/grad_J1b_ints.irp.f b/src/ao_many_one_e_ints/grad_J1b_ints.irp.f new file mode 100644 index 00000000..30e0acc8 --- /dev/null +++ b/src/ao_many_one_e_ints/grad_J1b_ints.irp.f @@ -0,0 +1,398 @@ + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb_size] + + implicit none + + List_all_comb_size = 2**nucl_num + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb, (nucl_num, List_all_comb_size)] + + implicit none + integer :: i, j + + if(nucl_num .gt. 32) then + print *, ' nucl_num = ', nucl_num, '> 32' + stop + endif + + List_all_comb = 0 + + do i = 0, List_all_comb_size-1 + do j = 0, nucl_num-1 + if (btest(i,j)) then + List_all_comb(j+1,i+1) = 1 + endif + enddo + enddo + +END_PROVIDER + +! --- + + BEGIN_PROVIDER [ double precision, List_all_j1b1s_coef, ( List_all_comb_size)] +&BEGIN_PROVIDER [ double precision, List_all_j1b1s_expo, ( List_all_comb_size)] +&BEGIN_PROVIDER [ double precision, List_all_j1b1s_cent, (3, List_all_comb_size)] + + implicit none + integer :: i, j, k, phase + double precision :: tmp_alphaj, tmp_alphak + + provide j1b_pen + + List_all_j1b1s_coef = 0.d0 + List_all_j1b1s_expo = 0.d0 + List_all_j1b1s_cent = 0.d0 + + do i = 1, List_all_comb_size + + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb(j,i)) * j1b_pen(j) + + List_all_j1b1s_expo(i) += tmp_alphaj + List_all_j1b1s_cent(1,i) += tmp_alphaj * nucl_coord(j,1) + List_all_j1b1s_cent(2,i) += tmp_alphaj * nucl_coord(j,2) + List_all_j1b1s_cent(3,i) += tmp_alphaj * nucl_coord(j,3) + + enddo + + ASSERT(List_all_j1b1s_expo(i) .gt. 0d0) + if(List_all_j1b1s_expo(i) .lt. 1d-10) cycle + + List_all_j1b1s_cent(1,i) = List_all_j1b1s_cent(1,i) / List_all_j1b1s_expo(i) + List_all_j1b1s_cent(2,i) = List_all_j1b1s_cent(2,i) / List_all_j1b1s_expo(i) + List_all_j1b1s_cent(3,i) = List_all_j1b1s_cent(3,i) / List_all_j1b1s_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_size + + do j = 2, nucl_num, 1 + tmp_alphaj = dble(List_all_comb(j,i)) * j1b_pen(j) + do k = 1, j-1, 1 + tmp_alphak = dble(List_all_comb(k,i)) * j1b_pen(k) + + List_all_j1b1s_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & + + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & + + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + enddo + enddo + + if(List_all_j1b1s_expo(i) .lt. 1d-10) cycle + + List_all_j1b1s_coef(i) = List_all_j1b1s_coef(i) / List_all_j1b1s_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_size + + phase = 0 + do j = 1, nucl_num + phase += List_all_comb(j,i) + enddo + + List_all_j1b1s_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_j1b1s_coef(i)) + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R| + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: r(3), int_mu, int_coulomb + double precision :: coef, beta, B_center(3) + double precision :: wall0, wall1 + double precision, allocatable :: tmp(:,:,:) + + double precision, external :: NAI_pol_mult_erf_ao_with1s + + provide mu_erf final_grid_points j1b_pen + call wall_time(wall0) + + v_ij_erf_rk_cst_mu_j1b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_size, final_grid_points, & + !$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, List_all_j1b1s_cent, & + !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) + + allocate( tmp(ao_num,ao_num,n_points_final_grid) ) + tmp = 0.d0 + + !$OMP DO + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i_1s = 1, List_all_comb_size + + coef = List_all_j1b1s_coef (i_1s) + beta = List_all_j1b1s_expo (i_1s) + B_center(1) = List_all_j1b1s_cent(1,i_1s) + B_center(2) = List_all_j1b1s_cent(2,i_1s) + B_center(3) = List_all_j1b1s_cent(3,i_1s) + + int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) + int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) + + tmp(j,i,ipoint) += coef * (int_mu - int_coulomb) + enddo + + enddo + enddo + enddo + !$OMP END DO + + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint) + enddo + enddo + enddo + !$OMP END CRITICAL + + deallocate( tmp ) + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_erf_rk_cst_mu_j1b', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)] + + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) + double precision :: wall0, wall1 + double precision, allocatable :: tmp(:,:,:,:) + + call wall_time(wall0) + + x_v_ij_erf_rk_cst_mu_tmp_j1b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_size, final_grid_points, & + !$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, List_all_j1b1s_cent, & + !$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b, mu_erf) + + allocate( tmp(3,ao_num,ao_num,n_points_final_grid) ) + tmp = 0.d0 + + !$OMP DO + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i_1s = 1, List_all_comb_size + + coef = List_all_j1b1s_coef (i_1s) + beta = List_all_j1b1s_expo (i_1s) + B_center(1) = List_all_j1b1s_cent(1,i_1s) + B_center(2) = List_all_j1b1s_cent(2,i_1s) + B_center(3) = List_all_j1b1s_cent(3,i_1s) + + 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) + + tmp(1,j,i,ipoint) += coef * (ints(1) - ints_coulomb(1)) + tmp(2,j,i,ipoint) += coef * (ints(2) - ints_coulomb(2)) + tmp(3,j,i,ipoint) += coef * (ints(3) - ints_coulomb(3)) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) += tmp(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) += tmp(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) += tmp(3,j,i,ipoint) + enddo + enddo + enddo + !$OMP END CRITICAL + + deallocate( tmp ) + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for x_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), int_fit, expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: wall0, wall1 + double precision, allocatable :: tmp(:,:,:) + + double precision, external :: overlap_gauss_r12_ao_with1s + + provide mu_erf final_grid_points j1b_pen + call wall_time(wall0) + + v_ij_u_cst_mu_j1b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$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_size, & + !$OMP final_grid_points, n_max_fit_slat, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, & + !$OMP List_all_j1b1s_cent, v_ij_u_cst_mu_j1b) + + allocate( tmp(ao_num,ao_num,n_points_final_grid) ) + tmp = 0.d0 + + !$OMP DO + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i_1s = 1, List_all_comb_size + + coef = List_all_j1b1s_coef (i_1s) + beta = List_all_j1b1s_expo (i_1s) + B_center(1) = List_all_j1b1s_cent(1,i_1s) + B_center(2) = List_all_j1b1s_cent(2,i_1s) + B_center(3) = List_all_j1b1s_cent(3,i_1s) + + do i_fit = 1, n_max_fit_slat + + expo_fit = expo_gauss_j_mu_x(i_fit) + coef_fit = coef_gauss_j_mu_x(i_fit) + int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + + tmp(j,i,ipoint) += coef * coef_fit * int_fit + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + v_ij_u_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint) + enddo + enddo + enddo + !$OMP END CRITICAL + + deallocate( tmp ) + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + v_ij_u_cst_mu_j1b(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0 + +END_PROVIDER + +! --- diff --git a/src/ao_many_one_e_ints/grad_related_ints.irp.f b/src/ao_many_one_e_ints/grad_related_ints.irp.f index c3c886f8..13fb1fc8 100644 --- a/src/ao_many_one_e_ints/grad_related_ints.irp.f +++ b/src/ao_many_one_e_ints/grad_related_ints.irp.f @@ -1,47 +1,64 @@ -BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, ( ao_num, ao_num,n_points_final_grid)] - implicit none - BEGIN_DOC -! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1)/|r - R| - END_DOC - integer :: i,j,ipoint - double precision :: mu,r(3),NAI_pol_mult_erf_ao - double precision :: int_mu, int_coulomb - provide mu_erf final_grid_points - double precision :: wall0, wall1 - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,mu,r,int_mu,int_coulomb) & - !$OMP SHARED (ao_num,n_points_final_grid,v_ij_erf_rk_cst_mu,final_grid_points,mu_erf) + +! --- + +BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1) / |r - R| + ! + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: r(3) + double precision :: int_mu, int_coulomb + double precision :: wall0, wall1 + + double precision :: NAI_pol_mult_erf_ao + + provide mu_erf final_grid_points + call wall_time(wall0) + + v_ij_erf_rk_cst_mu = 0.d0 + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint, r, int_mu, int_coulomb) & + !$OMP SHARED (ao_num, n_points_final_grid, v_ij_erf_rk_cst_mu, final_grid_points, mu_erf) !$OMP DO SCHEDULE (dynamic) - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - mu = mu_erf - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - int_mu = NAI_pol_mult_erf_ao(i,j,mu,r) - int_coulomb = NAI_pol_mult_erf_ao(i,j,1.d+9,r) - v_ij_erf_rk_cst_mu(j,i,ipoint)= (int_mu - int_coulomb ) - enddo + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + int_mu = NAI_pol_mult_erf_ao(i, j, mu_erf, r) + int_coulomb = NAI_pol_mult_erf_ao(i, j, 1.d+9, r) + + v_ij_erf_rk_cst_mu(j,i,ipoint) = int_mu - int_coulomb + enddo + enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, i-1 - v_ij_erf_rk_cst_mu(j,i,ipoint)= v_ij_erf_rk_cst_mu(i,j,ipoint) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + v_ij_erf_rk_cst_mu(j,i,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) + enddo enddo - enddo enddo + + call wall_time(wall1) + print*, 'wall time for v_ij_erf_rk_cst_mu ', wall1 - wall0 - call wall_time(wall1) - print*,'wall time for v_ij_erf_rk_cst_mu ',wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_grid, ao_num, ao_num)] implicit none BEGIN_DOC @@ -86,54 +103,62 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_gr print*,'wall time for v_ij_erf_rk_cst_mu_transp ',wall1 - wall0 END_PROVIDER +! --- + +BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint, m + double precision :: r(3), ints(3), ints_coulomb(3) + double precision :: wall0, wall1 + + call wall_time(wall0) -BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3,ao_num, ao_num,n_points_final_grid)] - implicit none - BEGIN_DOC -! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints(3),ints_coulomb(3) - double precision :: wall0, wall1 - call wall_time(wall0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,mu,r,ints,m,ints_coulomb) & + !$OMP PRIVATE (i,j,ipoint,r,ints,m,ints_coulomb) & !$OMP SHARED (ao_num,n_points_final_grid,x_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf) !$OMP DO SCHEDULE (dynamic) - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - mu = mu_erf - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - call NAI_pol_x_mult_erf_ao(i,j,mu,r,ints) - call NAI_pol_x_mult_erf_ao(i,j,1.d+9,r,ints_coulomb) - do m = 1, 3 - x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = ( ints(m) - ints_coulomb(m)) + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + call NAI_pol_x_mult_erf_ao(i, j, mu_erf, r, ints ) + call NAI_pol_x_mult_erf_ao(i, j, 1.d+9 , r, ints_coulomb) + + do m = 1, 3 + x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = (ints(m) - ints_coulomb(m)) + enddo + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, i-1 - do m = 1, 3 - x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint)= x_v_ij_erf_rk_cst_mu_tmp(m,i,j,ipoint) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + do m = 1, 3 + x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(m,i,j,ipoint) + enddo + enddo enddo - enddo enddo - enddo - call wall_time(wall1) - print*,'wall time for x_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0 + call wall_time(wall1) + print*,'wall time for x_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu, (ao_num, ao_num,n_points_final_grid,3)] implicit none BEGIN_DOC 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 42505194..d4ef4b28 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 @@ -1,3 +1,6 @@ + +! --- + subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) implicit none BEGIN_DOC @@ -15,142 +18,328 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) enddo end +! --- + +double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) -double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center) - implicit none BEGIN_DOC + ! ! Computes the following integral : - ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! END_DOC - integer, intent(in) :: i_ao,j_ao + + implicit none + integer, intent(in) :: i_ao, j_ao double precision, intent(in) :: mu_in, C_center(3) - integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in - double precision :: A_center(3), B_center(3),integral, alpha,beta + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in + double precision :: A_center(3), B_center(3), integral, alpha, beta + double precision :: NAI_pol_mult_erf - num_A = ao_nucl(i_ao) - power_A(1:3)= ao_power(i_ao,1:3) + + num_A = ao_nucl(i_ao) + power_A(1:3) = ao_power(i_ao,1:3) A_center(1:3) = nucl_coord(num_A,1:3) - num_B = ao_nucl(j_ao) - power_B(1:3)= ao_power(j_ao,1:3) + num_B = ao_nucl(j_ao) + power_B(1:3) = ao_power(j_ao,1:3) B_center(1:3) = nucl_coord(num_B,1:3) + n_pt_in = n_pt_max_integrals + NAI_pol_mult_erf_ao = 0.d0 do i = 1, ao_prim_num(i_ao) alpha = ao_expo_ordered_transp(i,i_ao) do j = 1, ao_prim_num(j_ao) beta = ao_expo_ordered_transp(j,j_ao) - integral = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in) - NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao) + + integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in) + + NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) enddo enddo -end +end function NAI_pol_mult_erf_ao +! --- +double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center) -double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in) BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + implicit none + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: beta, B_center(3) + double precision, intent(in) :: mu_in, C_center(3) + + integer :: i, j, power_A1(3), power_A2(3), n_pt_in + double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, integral + + double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao + + ASSERT(beta .lt. 0.d0) + if(beta .lt. 1d-10) then + NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) + return + endif + + power_A1(1:3) = ao_power(i_ao,1:3) + power_A2(1:3) = ao_power(j_ao,1:3) + + A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3) + A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3) + + n_pt_in = n_pt_max_integrals + + NAI_pol_mult_erf_ao_with1s = 0.d0 + do i = 1, ao_prim_num(i_ao) + alpha1 = ao_expo_ordered_transp(i,i_ao) + do j = 1, ao_prim_num(j_ao) + alpha2 = ao_expo_ordered_transp(j,j_ao) + + coef12 = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) + if(coef12 .lt. 1d-14) cycle + + integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & + , beta, B_center, C_center, n_pt_in, mu_in ) + + NAI_pol_mult_erf_ao_with1s += integral * coef12 + enddo + enddo + +end function NAI_pol_mult_erf_ao_with1s + +! --- + +double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + BEGIN_DOC + ! ! Computes the following integral : ! ! .. math:: ! ! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 ) - ! \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$. ! END_DOC - implicit none - integer, intent(in) :: n_pt_in - double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta,mu_in - integer, intent(in) :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt - double precision :: P_center(3) - - double precision :: d(0:n_pt_in),pouet,coeff,dist,const,pouet_2,factor - double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi - double precision :: V_e_n,const_factor,dist_integral,tmp - double precision :: accu,rint,p_inv,p,rho,p_inv_2 - integer :: n_pt_out,lmax include 'utils/constants.include.F' - p = alpha + beta - p_inv = 1.d0/p - p_inv_2 = 0.5d0 * p_inv - rho = alpha * beta * p_inv - dist = 0.d0 + implicit none + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A(3), power_B(3) + double precision, intent(in) :: C_center(3), A_center(3), B_center(3), alpha, beta, mu_in + + integer :: i, n_pt, n_pt_out + double precision :: P_center(3) + double precision :: d(0:n_pt_in), coeff, dist, const, factor + double precision :: const_factor, dist_integral + double precision :: accu, p_inv, p, rho, p_inv_2 + double precision :: p_new + + double precision :: rint + + p = alpha + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha * beta * p_inv + + dist = 0.d0 dist_integral = 0.d0 do i = 1, 3 - P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv - dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i)) - dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i)) + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i)) + dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) enddo - const_factor = dist*rho - if(const_factor > 80.d0)then + const_factor = dist * rho + if(const_factor > 80.d0) then NAI_pol_mult_erf = 0.d0 return endif - double precision :: p_new - p_new = mu_in/dsqrt(p+ mu_in * mu_in) - factor = dexp(-const_factor) - coeff = dtwo_pi * factor * p_inv * p_new - lmax = 20 - ! print*, "b" + p_new = mu_in / dsqrt(p + mu_in * mu_in) + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + + n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) ) + const = p * dist_integral * p_new * p_new + if(n_pt == 0) then + NAI_pol_mult_erf = coeff * rint(0, const) + return + endif + do i = 0, n_pt_in d(i) = 0.d0 enddo - n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) ) - const = p * dist_integral * p_new * p_new - if (n_pt == 0) then - pouet = rint(0,const) - NAI_pol_mult_erf = coeff * pouet - return - endif - ! call 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) p_new = p_new * p_new - call give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center) + call give_polynomial_mult_center_one_e_erf_opt( A_center, B_center, power_A, power_B, C_center & + , n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) - - if(n_pt_out<0)then + if(n_pt_out < 0) then NAI_pol_mult_erf = 0.d0 return endif - accu = 0.d0 ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i - do i =0 ,n_pt_out,2 - accu += d(i) * rint(i/2,const) + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) enddo NAI_pol_mult_erf = accu * coeff -end +end function NAI_pol_mult_erf + +! --- + +double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & + , beta, B_center, C_center, n_pt_in, mu_in ) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! .. math:: + ! + ! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2) + ! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2) + ! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2) + ! \exp(-\beta (r - B)^2) + ! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A1(3), power_A2(3) + double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3) + double precision, intent(in) :: alpha1, alpha2, beta, mu_in + + integer :: i, n_pt, n_pt_out + double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12 + double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor + double precision :: dist_integral + double precision :: d(0:n_pt_in), coeff, const, factor + double precision :: accu + double precision :: p_new + + double precision :: rint -subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,& - power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center) + ! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{K12} e^{-alpha12 (r - A12)^2} + alpha12 = alpha1 + alpha2 + alpha12_inv = 1.d0 / alpha12 + alpha12_inv_2 = 0.5d0 * alpha12_inv + rho12 = alpha1 * alpha2 * alpha12_inv + + dist12 = 0.d0 + do i = 1, 3 + A12_center(i) = (alpha1 * A1_center(i) + alpha2 * A2_center(i)) * alpha12_inv + dist12 += (A1_center(i) - A2_center(i)) * (A1_center(i) - A2_center(i)) + enddo + + const_factor12 = dist12 * rho12 + if(const_factor12 > 80.d0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + ! --- + + ! e^{K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{K} e^{-p (r - P)^2} + p = alpha12 + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha12 * beta * p_inv + + dist = 0.d0 + dist_integral = 0.d0 + do i = 1, 3 + P_center(i) = (alpha12 * A12_center(i) + beta * B_center(i)) * p_inv + dist += (A12_center(i) - B_center(i)) * (A12_center(i) - B_center(i)) + dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) + enddo + + const_factor = const_factor12 + dist * rho + if(const_factor > 80.d0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + ! --- + + p_new = mu_in / dsqrt(p + mu_in * mu_in) + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + + n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) ) + const = p * dist_integral * p_new * p_new + if(n_pt == 0) then + NAI_pol_mult_erf_with1s = coeff * rint(0, const) + return + endif + + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + p_new = p_new * p_new + + call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center & + , n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) + + if(n_pt_out < 0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) + enddo + NAI_pol_mult_erf_with1s = accu * coeff + +end function NAI_pol_mult_erf_with1s + +! --- + +subroutine give_polynomial_mult_center_one_e_erf_opt( A_center, B_center, power_A, power_B, C_center & + , n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) + 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 - double precision, intent(in) :: A_center(3), B_center(3),C_center(3),p,p_inv,p_inv_2,p_new,P_center(3) - double precision, intent(in) :: alpha,beta,mu_in - integer, intent(in) :: power_A(3), power_B(3) - integer :: a_x,b_x,a_y,b_y,a_z,b_z - double precision :: d(0:n_pt_in) - double precision :: d1(0:n_pt_in) - double precision :: d2(0:n_pt_in) - double precision :: d3(0:n_pt_in) - double precision :: accu + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A(3), power_B(3) + double precision, intent(in) :: A_center(3), B_center(3), C_center(3), p_inv_2, p_new, P_center(3) + integer, intent(out) :: n_pt_out + double precision, intent(out) :: d(0:n_pt_in) + + integer :: a_x, b_x, a_y, b_y, a_z, b_z + integer :: n_pt1, n_pt2, n_pt3, dim, i + integer :: n_pt_tmp + double precision :: d1(0:n_pt_in) + double precision :: d2(0:n_pt_in) + double precision :: d3(0:n_pt_in) + double precision :: accu + double precision :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2) + accu = 0.d0 ASSERT (n_pt_in > 1) - double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) R1x(0) = (P_center(1) - A_center(1)) R1x(1) = 0.d0 R1x(2) = -(P_center(1) - C_center(1))* p_new @@ -161,27 +350,22 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet !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(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 - d(i) = 0.d0 - enddo - do i = 0,n_pt_in + + do i = 0, n_pt_in + d (i) = 0.d0 d1(i) = 0.d0 - enddo - do i = 0,n_pt_in d2(i) = 0.d0 - enddo - do i = 0,n_pt_in d3(i) = 0.d0 enddo - integer :: n_pt1,n_pt2,n_pt3,dim,i + n_pt1 = n_pt_in n_pt2 = n_pt_in n_pt3 = n_pt_in a_x = power_A(1) b_x = power_B(1) - call I_x1_pol_mult_one_e(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in) + call I_x1_pol_mult_one_e(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in) if(n_pt1<0)then n_pt_out = -1 do i = 0,n_pt_in @@ -200,7 +384,7 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet !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) - call I_x1_pol_mult_one_e(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in) + call I_x1_pol_mult_one_e(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in) if(n_pt2<0)then n_pt_out = -1 do i = 0,n_pt_in @@ -209,41 +393,40 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet 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(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(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) - call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) - if(n_pt3<0)then + call I_x1_pol_mult_one_e(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in) + if(n_pt3 < 0) then n_pt_out = -1 do i = 0,n_pt_in d(i) = 0.d0 enddo return endif - integer :: n_pt_tmp + n_pt_tmp = 0 - call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp) - do i = 0,n_pt_tmp + call multiply_poly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp) + do i = 0, n_pt_tmp d1(i) = 0.d0 enddo n_pt_out = 0 - call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out) + call multiply_poly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out) do i = 0, n_pt_out d(i) = d1(i) enddo -end - +end subroutine give_polynomial_mult_center_one_e_erf_opt +! --- subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,& diff --git a/src/non_h_ints_mu/fit_j.irp.f b/src/non_h_ints_mu/fit_j.irp.f index 695ead7f..d359d489 100644 --- a/src/non_h_ints_mu/fit_j.irp.f +++ b/src/non_h_ints_mu/fit_j.irp.f @@ -12,33 +12,45 @@ BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ] 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)] - implicit none - BEGIN_DOC -! J(mu,r12) = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) is expressed as -! -! J(mu,r12) = 0.5/mu * F(r12*mu) where F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) -! -! F(x) is fitted by - 1/sqrt(pi) * exp(-alpha * x) exp(-beta*mu^2x^2) (see expo_j_xmu) -! -! The slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians -! -! See Appendix 2 of JCP 154, 084119 (2021) -! - END_DOC - integer :: i - double precision :: expos(n_max_fit_slat),alpha,beta - alpha = expo_j_xmu(1) * mu_erf - call expo_fit_slater_gam(alpha,expos) - beta = expo_j_xmu(2) * mu_erf**2.d0 - - do i = 1, n_max_fit_slat - expo_gauss_j_mu_x(i) = expos(i) + beta - coef_gauss_j_mu_x(i) = coef_fit_slat_gauss(i) / (2.d0 * mu_erf) * (- 1/dsqrt(dacos(-1.d0))) - enddo + + BEGIN_DOC + ! + ! J(mu,r12) = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) is expressed as + ! + ! J(mu,r12) = 0.5/mu * F(r12*mu) where F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) + ! + ! F(x) is fitted by - 1/sqrt(pi) * exp(-alpha * x) exp(-beta*mu^2x^2) (see expo_j_xmu) + ! + ! The slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians + ! + ! See Appendix 2 of JCP 154, 084119 (2021) + ! + END_DOC + + implicit none + integer :: i + double precision :: tmp + double precision :: expos(n_max_fit_slat), alpha, beta + + tmp = -0.5d0 / (mu_erf * sqrt(dacos(-1.d0))) + + alpha = expo_j_xmu(1) * mu_erf + call expo_fit_slater_gam(alpha, expos) + beta = expo_j_xmu(2) * mu_erf * mu_erf + + 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) + enddo + END_PROVIDER +! --- + double precision function F_x_j(x) implicit none BEGIN_DOC @@ -89,3 +101,6 @@ double precision function j_mu_fit_gauss(x) enddo end + +! --- + diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 068381b4..8f4883eb 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -1,30 +1,84 @@ -BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, ( ao_num, ao_num,n_points_final_grid,3)] - implicit none - BEGIN_DOC - ! grad_1_u_ij_mu(i,j,ipoint) = -1 * \int dr2 \grad_r1 u(r1,r2) \phi_i(r2) \phi_j(r2) - ! - ! where r1 = r(ipoint) - ! - ! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 (r1 - r2) (erf(mu * r12)-1)/2 r_12 \phi_i(r2) \phi_j(r2) - END_DOC - integer :: ipoint,i,j,m - double precision :: r(3) - do m = 1, 3 - do ipoint = 1, n_points_final_grid - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - do j = 1, ao_num - do i = 1, ao_num - grad_1_u_ij_mu(i,j,ipoint,m) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(m) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,m) + +! --- + +BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num,n_points_final_grid, 3)] + + BEGIN_DOC + ! + ! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) + ! = \int dr2 [(r1 - r2) (erf(mu * r12)-1)/2 r_12] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) + ! + ! where r1 = r(ipoint) + ! + END_DOC + + implicit none + integer :: ipoint, i, j, i_1s + double precision :: r(3) + double precision :: a, d, e, tmp1, tmp2, fact_r, fact_xyz(3) + + PROVIDE j1b_type j1b_pen + + if(j1b_type .eq. 3) then + + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + fact_r = 1.d0 + fact_xyz(1) = 1.d0 + fact_xyz(2) = 1.d0 + fact_xyz(3) = 1.d0 + do i_1s = 1, nucl_num + a = j1b_pen(i_1s) + d = (r(1) - nucl_coord(i_1s,1)) * (r(1) - nucl_coord(i_1s,1)) & + + (r(2) - nucl_coord(i_1s,2)) * (r(2) - nucl_coord(i_1s,2)) & + + (r(3) - nucl_coord(i_1s,3)) * (r(3) - nucl_coord(i_1s,3)) + e = dexp(-a*d) + + fact_r = fact_r * (1.d0 - e) + fact_xyz(1) = fact_xyz(1) * (2.d0 * a * (r(1) - nucl_coord(i_1s,1)) * e) + fact_xyz(2) = fact_xyz(2) * (2.d0 * a * (r(2) - nucl_coord(i_1s,2)) * e) + fact_xyz(3) = fact_xyz(3) * (2.d0 * a * (r(3) - nucl_coord(i_1s,3)) * e) + enddo + + do j = 1, ao_num + do i = 1, ao_num + + tmp1 = fact_r * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b (i,j,ipoint) + + grad_1_u_ij_mu(i,j,ipoint,1) = tmp1 * r(1) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + fact_xyz(1) * tmp2 + grad_1_u_ij_mu(i,j,ipoint,2) = tmp1 * r(2) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + fact_xyz(2) * tmp2 + grad_1_u_ij_mu(i,j,ipoint,3) = tmp1 * r(3) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + fact_xyz(3) * tmp2 + enddo + enddo enddo - enddo - enddo - enddo - grad_1_u_ij_mu *= 0.5d0 + + else + + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + grad_1_u_ij_mu(i,j,ipoint,1) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(1) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1) + grad_1_u_ij_mu(i,j,ipoint,2) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(2) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2) + grad_1_u_ij_mu(i,j,ipoint,3) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(3) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3) + enddo + enddo + enddo + + endif + + grad_1_u_ij_mu *= 0.5d0 END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)] implicit none BEGIN_DOC