diff --git a/external/irpf90 b/external/irpf90 index 43160c60..4ab1b175 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b +Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 diff --git a/plugins/local/extra_basis_int/coul_1s.irp.f b/plugins/local/extra_basis_int/coul_1s.irp.f index 964ed36a..62ec4331 100644 --- a/plugins/local/extra_basis_int/coul_1s.irp.f +++ b/plugins/local/extra_basis_int/coul_1s.irp.f @@ -11,6 +11,11 @@ double precision function coul_full_ao_pq_r_1s(p,q,R,R_p,R_q) double precision, intent(in) :: R(3),R_p(3),R_q(3) integer, intent(in) :: p,q double precision :: coef,dist,P_pq(3),coefaos + if(.not.ao_extra_only_1s)then + print*,'You are using a function assuming that the extra basis is fitted on 1s functions' + print*,'But this is not the case apparently ... stopping' + stop + endif coefaos= ao_extra_coef_normalized(p,1) * ao_extra_coef_normalized(q,1) coef = inv_pi_gamma_pq_3_2_ao_extra(p,q) * E_pq_ao_extra(p,q) P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q @@ -40,6 +45,11 @@ double precision function coul_pq_r_1s(p,q,R,R_p,R_q) double precision, intent(in) :: R(3),R_p(3),R_q(3) integer, intent(in) :: p,q double precision :: dist,P_pq(3) + if(.not.ao_extra_only_1s)then + print*,'You are using a function assuming that the extra basis is fitted on 1s functions' + print*,'But this is not the case apparently ... stopping' + stop + endif P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q P_pq = P_pq * inv_gamma_pq_ao_extra(q,p) dist = (P_pq(1)-R(1)) * (P_pq(1)-R(1)) diff --git a/plugins/local/extra_basis_int/pot_ao_ints.irp.f b/plugins/local/extra_basis_int/pot_ao_ints.irp.f index 1828f980..5a7d7580 100644 --- a/plugins/local/extra_basis_int/pot_ao_ints.irp.f +++ b/plugins/local/extra_basis_int/pot_ao_ints.irp.f @@ -3,23 +3,45 @@ BEGIN_PROVIDER [ double precision, pot_vne_A_extra_basis, (ao_extra_num,ao_extra BEGIN_DOC ! ! Computes the following integral : - ! $\sum_{R in nuclei} -Z $ + ! $\sum_{R in the USUAL nuclei} -Z $ ! ! - ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the extra basis + ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis END_DOC - integer :: mu,nu,k_nucl - double precision :: mu_in, R_nucl(3),charge_nucl, integral - double precision :: NAI_pol_mult_erf_ao_extra - mu_in = 10.d0**10 + integer :: mu,nu + double precision :: v_nucl_extra_ao pot_vne_A_extra_basis = 0.d0 do mu = 1, ao_extra_num do nu = 1, ao_extra_num - do k_nucl = 1, nucl_num - R_nucl(1:3) = nucl_coord_transp(1:3,k_nucl) - charge_nucl = nucl_charge(k_nucl) - integral = NAI_pol_mult_erf_ao_extra(mu, nu, mu_in, R_nucl) - pot_vne_A_extra_basis(nu,mu) += -integral * charge_nucl + pot_vne_A_extra_basis(nu,mu)= v_nucl_extra_ao(mu,nu) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, pot_vne_extra_basis, (ao_num,ao_num)] + implicit none + BEGIN_DOC + ! + ! Computes the following integral : + ! $\sum_{R in EXTRA nuclei} -Z $ + ! + ! + ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the USUAL basis + END_DOC + integer :: mu,nu,k_nucl + double precision :: mu_in, R_nucl(3),charge_nucl, integral + double precision :: NAI_pol_mult_erf_ao + mu_in = 10.d0**10 + pot_vne_extra_basis = 0.d0 + do mu = 1, ao_num + do nu = 1, ao_num + do k_nucl = 1, extra_nucl_num + R_nucl(1:3) = extra_nucl_coord_transp(1:3,k_nucl) + charge_nucl = extra_nucl_charge(k_nucl) + integral = NAI_pol_mult_erf_ao(mu, nu, mu_in, R_nucl) + pot_vne_extra_basis(nu,mu) += -integral * charge_nucl enddo enddo enddo diff --git a/plugins/local/extra_basis_int/v_mixed_extra.irp.f b/plugins/local/extra_basis_int/v_mixed_extra.irp.f index ac856692..15958ab5 100644 --- a/plugins/local/extra_basis_int/v_mixed_extra.irp.f +++ b/plugins/local/extra_basis_int/v_mixed_extra.irp.f @@ -1,3 +1,6 @@ +!!! TODO:: optimize when "ao_extra_only_1s" is True + + double precision function v_extra_nucl_extra_ao(i_ao,j_ao) implicit none BEGIN_DOC @@ -23,6 +26,30 @@ double precision function v_extra_nucl_extra_ao(i_ao,j_ao) enddo end +double precision function v_extra_nucl_ao(i_ao,j_ao) + implicit none + BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne(r)$. + ! + ! + ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the REGULAR basis + ! + ! and v_ne(r) is the Coulomb potential coming from the EXTRA nuclei + END_DOC + integer, intent(in) ::i_ao,j_ao + integer :: i + double precision :: mu_in, coord(3),charge, integral + double precision :: NAI_pol_mult_erf_ao + mu_in = 1.d+10 + do i = 1, extra_nucl_num + coord(1:3) = extra_nucl_coord_transp(1:3,i) + charge = extra_nucl_charge(i) + v_extra_nucl_ao += -NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, coord) * charge + enddo +end + double precision function v_nucl_extra_ao(i_ao,j_ao) implicit none diff --git a/src/ao_extra_basis/fit_1s_basis.irp.f b/src/ao_extra_basis/fit_1s_basis.irp.f index ef09d5b2..9c788173 100644 --- a/src/ao_extra_basis/fit_1s_basis.irp.f +++ b/src/ao_extra_basis/fit_1s_basis.irp.f @@ -31,6 +31,7 @@ program fit_1s_basis call ezfio_set_extra_nuclei_extra_nucl_label(new_nucl_label_1s) ! call ezfio_set_ao_extra_basis_ao_extra_num(n_func_tot) + call ezfio_set_ao_extra_basis_ao_extra_only_1s(.True.) call ezfio_set_ao_extra_basis_ao_extra_center(ao_extra_center) call ezfio_set_ao_extra_basis_ao_extra_nucl(new_ao_nucl_1s) call ezfio_set_ao_extra_basis_ao_extra_prim_num(new_ao_prim_num_1s)