diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 7ee827d3..8ce54f3d 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -846,6 +846,71 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_cart_num,ao_sphe_num)] + implicit none + BEGIN_DOC +! Coefficients to go from cartesian to spherical coordinates in the current +! basis set +! +! S_cart^-1 + END_DOC + integer :: i + integer, external :: ao_power_index + integer :: ibegin,j,k + integer :: prev, ao_sphe_count + prev = 0 + ao_cart_to_sphe_coef(:,:) = 0.d0 + + if (ao_cartesian) then + ! Identity matrix + do i=1,ao_sphe_num + ao_cart_to_sphe_coef(i,i) = 1.d0 + enddo + + else + ! Assume order provided by ao_power_index + i = 1 + ao_sphe_count = 0 + do while (i <= ao_num) + select case ( ao_l(i) ) + case (0) + ao_sphe_count += 1 + ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 + i += 1 + BEGIN_TEMPLATE + case ($SHELL) + if (ao_power(i,1) == $SHELL) then + do k=1,size(cart_to_sphe_$SHELL,2) + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) + enddo + enddo + i += size(cart_to_sphe_$SHELL,1) + ao_sphe_count += size(cart_to_sphe_$SHELL,2) + endif + SUBST [ SHELL ] + 1;; + 2;; + 3;; + 4;; + 5;; + 6;; + 7;; + 8;; + 9;; + END_TEMPLATE + case default + stop 'Error in ao_cart_to_sphe : angular momentum too high' + end select + enddo + + endif + + if (ao_sphe_count /= ao_sphe_num) then + call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") + endif +END_PROVIDER + diff --git a/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f index 6919b51f..e2368278 100644 --- a/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ b/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -26,7 +26,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_ - ao_cart_integrals_n_e_cgtos = 0.d0 + ao_cart_coul_n_e_cgtos = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & @@ -37,7 +37,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_ !$OMP ao_cart_power, nucl_num, nucl_charge, n_pt_max_integrals, & !$OMP ao_cart_expo_cgtos_ord_transp, ao_cart_coef_cgtos_norm_ord_transp, & !$OMP ao_cart_expo_pw_ord_transp, ao_cart_expo_phase_ord_transp, & - !$OMP ao_cart_integrals_n_e_cgtos) + !$OMP ao_cart_coul_n_e_cgtos) !$OMP DO SCHEDULE (dynamic) do j = 1, ao_cart_num @@ -91,7 +91,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_ c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2) enddo - ao_cart_integrals_n_e_cgtos(i,j) += c * ao_cart_coef_cgtos_norm_ord_transp(n,j) & + ao_cart_coul_n_e_cgtos(i,j) += c * ao_cart_coef_cgtos_norm_ord_transp(n,j) & * ao_cart_coef_cgtos_norm_ord_transp(l,i) enddo enddo diff --git a/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f b/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f index a81206c1..93fec121 100644 --- a/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f +++ b/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f @@ -1,7 +1,7 @@ ! --- -subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) +subroutine give_all_erf_kl_ao_cart(integrals_ao,mu_in,C_center) implicit none BEGIN_DOC ! Subroutine that returns all integrals over $r$ of type @@ -9,18 +9,18 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) END_DOC double precision, intent(in) :: mu_in,C_center(3) double precision, intent(out) :: integrals_ao(ao_cart_num,ao_cart_num) - double precision :: NAI_pol_mult_erf_ao + double precision :: NAI_pol_mult_erf_ao_cart integer :: i,j,l,k,m do k = 1, ao_cart_num do m = 1, ao_cart_num - integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center) + integrals_ao(m,k) = NAI_pol_mult_erf_ao_cart(m,k,mu_in,C_center) enddo 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_cart(i_ao, j_ao, mu_in, C_center) BEGIN_DOC ! diff --git a/src/ao_cart_one_e_ints/pot_ao_ints.irp.f b/src/ao_cart_one_e_ints/pot_ao_ints.irp.f index 4c59afb7..7c79e0b2 100644 --- a/src/ao_cart_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_cart_one_e_ints/pot_ao_ints.irp.f @@ -1,7 +1,7 @@ ! --- -BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_num)] +BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e, (ao_cart_num,ao_cart_num)] BEGIN_DOC ! Nucleus-electron interaction, in the |AO| basis set. @@ -18,11 +18,11 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n double precision :: A_center(3),B_center(3),C_center(3) double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - ao_cart_integrals_n_e = 0.d0 + ao_cart_coul_n_e = 0.d0 - if (read_ao_cart_integrals_n_e) then + if (read_ao_cart_coul_n_e) then - call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e) + call ezfio_get_ao_cart_one_e_ints_ao_cart_coul_n_e(ao_cart_coul_n_e) print *, 'AO N-e integrals read from disk' else @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n do j = 1, ao_cart_num do i = 1, ao_cart_num - ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e_cgtos(i,j) + ao_cart_coul_n_e(i,j) = ao_cart_coul_n_e_cgtos(i,j) enddo enddo @@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& !$OMP num_A,num_B,Z,c,c1,n_pt_in) & !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& - !$OMP n_pt_max_integrals,ao_cart_integrals_n_e,nucl_num,nucl_charge) + !$OMP n_pt_max_integrals,ao_cart_coul_n_e,nucl_num,nucl_charge) n_pt_in = n_pt_max_integrals @@ -86,7 +86,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n c = c - Z * c1 enddo - ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e(i,j) & + ao_cart_coul_n_e(i,j) = ao_cart_coul_n_e(i,j) & + ao_cart_coef_normalized_ordered_transp(l,j) & * ao_cart_coef_normalized_ordered_transp(m,i) * c enddo @@ -99,25 +99,17 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n endif - - IF(do_pseudo) THEN - ao_cart_integrals_n_e += ao_cart_pseudo_integrals - ENDIF - IF(point_charges) THEN - ao_cart_integrals_n_e += ao_cart_integrals_pt_chrg - ENDIF - endif - if (write_ao_cart_integrals_n_e) then - call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e) + if (write_ao_cart_coul_n_e) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_coul_n_e(ao_cart_coul_n_e) print *, 'AO N-e integrals written to disk' endif END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_cart_num)] +BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e_imag, (ao_cart_num,ao_cart_num)] BEGIN_DOC ! Nucleus-electron interaction, in the |AO| basis set. ! @@ -131,8 +123,8 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_c integer :: i,j,k,l,n_pt_in,m double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - if (read_ao_cart_integrals_n_e) then - call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e_imag(ao_cart_integrals_n_e_imag) + if (read_ao_cart_coul_n_e) then + call ezfio_get_ao_cart_one_e_ints_ao_cart_coul_n_e_imag(ao_cart_coul_n_e_imag) print *, 'AO N-e integrals read from disk' else print *, irp_here, ': Not yet implemented' @@ -140,7 +132,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_c END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num,ao_cart_num,nucl_num)] +BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e_per_atom, (ao_cart_num,ao_cart_num,nucl_num)] BEGIN_DOC ! Nucleus-electron interaction in the |AO| basis set, per atom A. ! @@ -154,14 +146,14 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num, integer :: i,j,k,l,n_pt_in,m double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - ao_cart_integrals_n_e_per_atom = 0.d0 + ao_cart_coul_n_e_per_atom = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,power_A,power_B,& !$OMP num_A,num_B,c,n_pt_in,C_center) & !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& - !$OMP n_pt_max_integrals,ao_cart_integrals_n_e_per_atom,nucl_num) + !$OMP n_pt_max_integrals,ao_cart_coul_n_e_per_atom,nucl_num) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (dynamic) @@ -197,7 +189,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num, * ao_cart_coef_normalized_ordered_transp(m,i) enddo enddo - ao_cart_integrals_n_e_per_atom(i,j,k) = -c + ao_cart_coul_n_e_per_atom(i,j,k) = -c enddo enddo enddo diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 8d4fff57..74bb5988 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -16,7 +16,6 @@ doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write interface: ezfio,provider,ocaml default: None - [ao_integrals_kinetic] type: double precision doc: Kinetic energy integrals in |AO| basis set diff --git a/src/ao_one_e_ints/NEED b/src/ao_one_e_ints/NEED index 61d23b1e..b7d0ccdb 100644 --- a/src/ao_one_e_ints/NEED +++ b/src/ao_one_e_ints/NEED @@ -1,2 +1,3 @@ ao_basis +ao_cart_one_e_ints pseudo diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f index 8b0edfbc..2c6b8e7e 100644 --- a/src/ao_one_e_ints/ao_one_e_ints.irp.f +++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f @@ -45,19 +45,3 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals,(ao_sphe_num,ao_sphe_num)] -&BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_diag,(ao_sphe_num)] - implicit none - integer :: i,j,n,l - BEGIN_DOC - ! One-electron Hamiltonian in the spherical |AO| basis. - END_DOC - - ao_sphe_one_e_integrals = ao_sphe_integrals_n_e + ao_sphe_kinetic_integrals - - do j = 1, ao_num - ao_sphe_one_e_integrals_diag(j) = ao_sphe_one_e_integrals(j,j) - enddo - -END_PROVIDER - diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index e015c89e..48268c73 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -1,74 +1,3 @@ - BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)] -&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)] - implicit none - BEGIN_DOC -! Coefficients to go from cartesian to spherical coordinates in the current -! basis set -! -! S_cart^-1 - END_DOC - integer :: i - integer, external :: ao_power_index - integer :: ibegin,j,k - integer :: prev, ao_sphe_count - prev = 0 - ao_cart_to_sphe_coef(:,:) = 0.d0 - ao_cart_to_sphe_normalization(:) = 1.d0 - - if (ao_cartesian) then - ! Identity matrix - do i=1,ao_sphe_num - ao_cart_to_sphe_coef(i,i) = 1.d0 - enddo - - else - ! Assume order provided by ao_power_index - i = 1 - ao_sphe_count = 0 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_sphe_count += 1 - ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 - ao_cart_to_sphe_normalization(i) = 1.d0 - i += 1 - BEGIN_TEMPLATE - case ($SHELL) - if (ao_power(i,1) == $SHELL) then - do k=1,size(cart_to_sphe_$SHELL,2) - do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) - enddo - enddo - do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j) - enddo - i += size(cart_to_sphe_$SHELL,1) - ao_sphe_count += size(cart_to_sphe_$SHELL,2) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in ao_cart_to_sphe : angular momentum too high' - end select - enddo - - endif - - if (ao_sphe_count /= ao_sphe_num) then - call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") - endif -END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] implicit none BEGIN_DOC @@ -146,63 +75,10 @@ END_PROVIDER enddo call write_double(6, lin_dep_cutoff, "Linear dependencies cut-off") - if (ao_cartesian) then - - ao_ortho_canonical_num = ao_num - call ortho_canonical(ao_overlap,size(ao_overlap,1), & - ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & - ao_ortho_canonical_num, lin_dep_cutoff) - - - else - - double precision, allocatable :: S(:,:) - - allocate(S(ao_sphe_num,ao_sphe_num)) - S = 0.d0 - do i=1,ao_sphe_num - S(i,i) = 1.d0 - enddo - - ao_ortho_canonical_num = ao_sphe_num - call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & - ao_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff) - - call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_sphe_num, 1.d0, & - ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & - S, size(S,1), & - 0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1)) - - deallocate(S) - endif - - FREE ao_overlap + ao_ortho_canonical_num = ao_num + call ortho_canonical(ao_overlap,size(ao_overlap,1), & + ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & + ao_ortho_canonical_num, lin_dep_cutoff) END_PROVIDER -BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] - implicit none - BEGIN_DOC -! overlap matrix of the ao_ortho_canonical. -! Expected to be the Identity - END_DOC - integer :: i,j,k,l - double precision :: c - do j=1, ao_ortho_canonical_num - do i=1, ao_ortho_canonical_num - ao_ortho_canonical_overlap(i,j) = 0.d0 - enddo - enddo - do j=1, ao_ortho_canonical_num - do k=1, ao_num - c = 0.d0 - do l=1, ao_num - c += ao_ortho_canonical_coef(l,j) * ao_overlap(l,k) - enddo - do i=1, ao_ortho_canonical_num - ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c - enddo - enddo - enddo -END_PROVIDER - diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 17eb7cbc..00cc6718 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -2,9 +2,6 @@ ! --- BEGIN_PROVIDER [double precision, ao_overlap , (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_x, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_y, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_z, (ao_num, ao_num)] BEGIN_DOC ! Overlap between atomic basis functions: @@ -19,78 +16,12 @@ double precision :: A_center(3), B_center(3) ao_overlap = 0.d0 - ao_overlap_x = 0.d0 - ao_overlap_y = 0.d0 - ao_overlap_z = 0.d0 if(read_ao_integrals_overlap) then - - call ezfio_get_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num)) - print *, 'AO overlap integrals read from disk' - + call ezfio_get_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num)) + print *, 'AO overlap integrals read from disk' else - - if(use_cgtos) then - - do j = 1, ao_num - do i = 1, ao_num - ao_overlap (i,j) = ao_overlap_cgtos (i,j) - ao_overlap_x(i,j) = ao_overlap_cgtos_x(i,j) - ao_overlap_y(i,j) = ao_overlap_cgtos_y(i,j) - ao_overlap_z(i,j) = ao_overlap_cgtos_z(i,j) - enddo - enddo - - else - - dim1=100 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_x,overlap_y, overlap_z, overlap, & - !$OMP alpha, beta,i,j,n,l,c) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) - ao_overlap(i,j) += c * overlap - if(isnan(ao_overlap(i,j)))then - print*,'i,j',i,j - print*,'l,n',l,n - print*,'c,overlap',c,overlap - print*,overlap_x,overlap_y,overlap_z - stop - endif - ao_overlap_x(i,j) += c * overlap_x - ao_overlap_y(i,j) += c * overlap_y - ao_overlap_z(i,j) += c * overlap_z - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - - endif - + call ao_cart_to_ao_basis(ao_cart_overlap, ao_cart_num, ao_overlap, ao_num) endif if (write_ao_integrals_overlap) then @@ -151,45 +82,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs, (ao_num, ao_num) ] enddo else - - dim1=100 - lower_exp_val = 40.d0 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B, & - !$OMP overlap_x,overlap_y, overlap_z, & - !$OMP alpha, beta,i,j,dx) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,& - !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - ao_overlap_abs(i,j)= 0.d0 - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) - call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1) - call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1) - ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO + print*,'todo ! numerical integration on DFT grid !' + stop endif @@ -308,26 +202,3 @@ BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ] END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! |AO| overlap matrix in the spherical basis set - END_DOC - double precision, allocatable :: tmp(:,:) - allocate (tmp(ao_sphe_num,ao_num)) - - call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & - ao_overlap,size(ao_overlap,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - tmp, size(tmp,1), & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & - ao_sphe_overlap,size(ao_sphe_overlap,1)) - - deallocate(tmp) - -END_PROVIDER - diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f deleted file mode 100644 index 83c782b7..00000000 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ /dev/null @@ -1,124 +0,0 @@ -! --- - - BEGIN_PROVIDER [double precision, ao_overlap_cgtos, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_x, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_y, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_z, (ao_num, ao_num)] - - implicit none - - integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) - double precision :: c, overlap, overlap_x, overlap_y, overlap_z - double precision :: KA2(3), phiA(3) - double precision :: KB2(3), phiB(3) - complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) - complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) - complex*16 :: C1(1:4), C2(1:4) - complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1 - complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2 - - ao_overlap_cgtos = 0.d0 - ao_overlap_cgtos_x = 0.d0 - ao_overlap_cgtos_y = 0.d0 - ao_overlap_cgtos_z = 0.d0 - - dim1 = 100 - - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & - !$OMP alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, & - !$OMP beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB, & - !$OMP overlap_x , overlap_y , overlap_z , overlap, & - !$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, & - !$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) & - !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & - !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, & - !$OMP ao_overlap_cgtos) - - do j = 1, ao_num - - jj = ao_nucl(j) - power_A(1) = ao_power(j,1) - power_A(2) = ao_power(j,2) - power_A(3) = ao_power(j,3) - - do i = 1, ao_num - - ii = ao_nucl(i) - power_B(1) = ao_power(i,1) - power_B(2) = ao_power(i,2) - power_B(3) = ao_power(i,3) - - do n = 1, ao_prim_num(j) - - alpha = ao_expo_cgtos_ord_transp(n,j) - alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 - phiA(m) = ao_expo_phase_ord_transp(m,n,j) - Ap_center(m) = nucl_coord(jj,m) - Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) - KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j) - enddo - - do l = 1, ao_prim_num(i) - - beta = ao_expo_cgtos_ord_transp(l,i) - beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 - phiB(m) = ao_expo_phase_ord_transp(m,l,i) - Bp_center(m) = nucl_coord(ii,m) - Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) - KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i) - enddo - - c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i) - - C1(1) = zexp((0.d0, 1.d0) * (-phiA(1) - phiB(1)) - 0.25d0 * (alpha_inv * KA2(1) + beta_inv * KB2(1))) - C1(2) = zexp((0.d0, 1.d0) * (-phiA(2) - phiB(2)) - 0.25d0 * (alpha_inv * KA2(2) + beta_inv * KB2(2))) - C1(3) = zexp((0.d0, 1.d0) * (-phiA(3) - phiB(3)) - 0.25d0 * (alpha_inv * KA2(3) + beta_inv * KB2(3))) - C1(4) = C1(1) * C1(2) * C1(3) - - C2(1) = zexp((0.d0, 1.d0) * (phiA(1) - phiB(1)) - 0.25d0 * (conjg(alpha_inv) * KA2(1) + beta_inv * KB2(1))) - C2(2) = zexp((0.d0, 1.d0) * (phiA(2) - phiB(2)) - 0.25d0 * (conjg(alpha_inv) * KA2(2) + beta_inv * KB2(2))) - C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3))) - C2(4) = C2(1) * C2(2) * C2(3) - - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) - - call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, & - conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) - - overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2) - overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2) - overlap_z = 2.d0 * real(C1(3) * overlap_z1 + C2(3) * overlap_z2) - overlap = 2.d0 * real(C1(4) * overlap1 + C2(4) * overlap2 ) - - ao_overlap_cgtos(i,j) = ao_overlap_cgtos(i,j) + c * overlap - - if(isnan(ao_overlap_cgtos(i,j))) then - print*,'i, j', i, j - print*,'l, n', l, n - print*,'c, overlap', c, overlap - print*, overlap_x, overlap_y, overlap_z - stop - endif - - ao_overlap_cgtos_x(i,j) = ao_overlap_cgtos_x(i,j) + c * overlap_x - ao_overlap_cgtos_y(i,j) = ao_overlap_cgtos_y(i,j) + c * overlap_y - ao_overlap_cgtos_z(i,j) = ao_overlap_cgtos_z(i,j) + c * overlap_z - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - -! --- - - - diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f index 1e3e684d..67b0e04f 100644 --- a/src/ao_one_e_ints/kin_ao_ints.irp.f +++ b/src/ao_one_e_ints/kin_ao_ints.irp.f @@ -16,123 +16,11 @@ END_DOC implicit none - integer :: i, j, n, l, dim1, power_A(3), power_B(3) - double precision :: overlap, overlap_y, overlap_z - double precision :: overlap_x0, overlap_y0, overlap_z0 - double precision :: alpha, beta, c - double precision :: A_center(3), B_center(3) - double precision :: d_a_2,d_2 + integer :: i, j - if(use_cgtos) then - - do j = 1, ao_num - do i = 1, ao_num - ao_deriv2_x(i,j) = ao_deriv2_cgtos_x(i,j) - ao_deriv2_y(i,j) = ao_deriv2_cgtos_y(i,j) - ao_deriv2_z(i,j) = ao_deriv2_cgtos_z(i,j) - enddo - enddo - - else - - dim1=100 - - ! -- Dummy call to provide everything - A_center(:) = 0.d0 - B_center(:) = 1.d0 - alpha = 1.d0 - beta = .1d0 - power_A = 1 - power_B = 0 - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1) - ! -- - - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_y, overlap_z, overlap, & - !$OMP alpha, beta, n, l, i,j,c,d_a_2,d_2,deriv_tmp, & - !$OMP overlap_x0,overlap_y0,overlap_z0) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - ao_deriv2_x(i,j)= 0.d0 - ao_deriv2_y(i,j)= 0.d0 - ao_deriv2_z(i,j)= 0.d0 - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1) - c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) - - power_A(1) = power_A(1)-2 - if (power_A(1)>-1) then - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_a_2,overlap_y,overlap_z,overlap,dim1) - else - d_a_2 = 0.d0 - endif - power_A(1) = power_A(1)+4 - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_2,overlap_y,overlap_z,overlap,dim1) - power_A(1) = power_A(1)-2 - - double precision :: deriv_tmp - deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(1)) +1.d0) * overlap_x0 & - +dble(power_A(1)) * (dble(power_A(1))-1.d0) * d_a_2 & - +4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0 - - ao_deriv2_x(i,j) += c*deriv_tmp - power_A(2) = power_A(2)-2 - if (power_A(2)>-1) then - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1) - else - d_a_2 = 0.d0 - endif - power_A(2) = power_A(2)+4 - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1) - power_A(2) = power_A(2)-2 - - deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(2)) +1.d0 ) * overlap_y0 & - +dble(power_A(2)) * (dble(power_A(2))-1.d0) * d_a_2 & - +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0 - ao_deriv2_y(i,j) += c*deriv_tmp - - power_A(3) = power_A(3)-2 - if (power_A(3)>-1) then - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1) - else - d_a_2 = 0.d0 - endif - power_A(3) = power_A(3)+4 - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1) - power_A(3) = power_A(3)-2 - - deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(3)) +1.d0 ) * overlap_z0 & - +dble(power_A(3)) * (dble(power_A(3))-1.d0) * d_a_2 & - +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0 - ao_deriv2_z(i,j) += c*deriv_tmp - - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - - endif + call ao_cart_to_ao_basis(ao_cart_deriv2_x, ao_cart_num, ao_deriv2_x, ao_num) + call ao_cart_to_ao_basis(ao_cart_deriv2_y, ao_cart_num, ao_deriv2_y, ao_num) + call ao_cart_to_ao_basis(ao_cart_deriv2_z, ao_cart_num, ao_deriv2_z, ao_num) END_PROVIDER @@ -190,25 +78,3 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)] endif END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_sphe_kinetic_integrals, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! |AO| kinetic inntegrals matrix in the spherical basis set - END_DOC - double precision, allocatable :: tmp(:,:) - allocate (tmp(ao_sphe_num,ao_num)) - - call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & - ao_kinetic_integrals,size(ao_kinetic_integrals,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - tmp, size(tmp,1), & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & - ao_sphe_kinetic_integrals,size(ao_sphe_kinetic_integrals,1)) - - deallocate(tmp) - -END_PROVIDER diff --git a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f deleted file mode 100644 index 7240f4a5..00000000 --- a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ /dev/null @@ -1,558 +0,0 @@ - -! --- - -BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] - - BEGIN_DOC - ! - ! Nucleus-electron interaction, in the cgtos |AO| basis set. - ! - ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` - ! - END_DOC - - implicit none - - integer :: power_A(3), power_B(3) - integer :: i, j, k, l, m, n, ii, jj - double precision :: c, Z, C_center(3) - double precision :: phiA, KA2 - double precision :: phiB, KB2 - complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) - complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) - complex*16 :: C1, C2, I1, I2 - - complex*16, external :: NAI_pol_mult_cgtos - - - - ao_integrals_n_e_cgtos = 0.d0 - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, C1, C2, I1, I2, & - !$OMP alpha, alpha_inv, Ae_center, Ap_center, phiA, KA2, power_A, & - !$OMP beta, beta_inv, Be_center, Bp_center, phiB, KB2, power_B) & - !$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, & - !$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, & - !$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_integrals_n_e_cgtos) - !$OMP DO SCHEDULE (dynamic) - - do j = 1, ao_num - - jj = ao_nucl(j) - power_A(1:3) = ao_power(j,1:3) - - do i = 1, ao_num - - ii = ao_nucl(i) - power_B(1:3) = ao_power(i,1:3) - - do n = 1, ao_prim_num(j) - - alpha = ao_expo_cgtos_ord_transp(n,j) - alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 - Ap_center(m) = nucl_coord(jj,m) - Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) - enddo - phiA = ao_expo_phase_ord_transp(4,n,j) - KA2 = ao_expo_pw_ord_transp(4,n,j) - - do l = 1, ao_prim_num(i) - - beta = ao_expo_cgtos_ord_transp(l,i) - beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 - Bp_center(m) = nucl_coord(ii,m) - Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) - enddo - phiB = ao_expo_phase_ord_transp(4,l,i) - KB2 = ao_expo_pw_ord_transp(4,l,i) - - C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) - C2 = zexp((0.d0, 1.d0) * ( phiA - phiB) - 0.25d0 * (conjg(alpha_inv) * KA2 + beta_inv * KB2)) - - c = 0.d0 - do k = 1, nucl_num - - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - I1 = NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & - Ap_center, Bp_center, C_center, n_pt_max_integrals) - - I2 = NAI_pol_mult_cgtos(conjg(Ae_center), Be_center, power_A, power_B, conjg(alpha), beta, & - conjg(Ap_center), Bp_center, C_center, n_pt_max_integrals) - - c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2) - enddo - - ao_integrals_n_e_cgtos(i,j) += c * ao_coef_cgtos_norm_ord_transp(n,j) & - * ao_coef_cgtos_norm_ord_transp(l,i) - enddo - enddo - enddo - enddo - - !$OMP END DO - !$OMP END PARALLEL - -END_PROVIDER - -! --- - -complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & - Ap_center, Bp_center, C_center, n_pt_in) - - BEGIN_DOC - ! - ! Computes the electron-nucleus attraction with two primitves cgtos. - ! - ! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle` - ! - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, power_A(3), power_B(3) - double precision, intent(in) :: C_center(3) - complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3) - complex*16, intent(in) :: beta, Be_center(3), Bp_center(3) - - integer :: i, n_pt, n_pt_out - double precision :: dist_AB, dist_AC - complex*16 :: p, p_inv, rho, dist, dist_integral, const, const_factor, coeff, factor - complex*16 :: P_center(3) - complex*16 :: d(0:n_pt_in) - - complex*16, external :: V_n_e_cgtos - complex*16, external :: crint_sum - complex*16, external :: crint_1 - - - - dist_AB = 0.d0 - dist_AC = 0.d0 - do i = 1, 3 - dist_AB += abs(Ae_center(i) - Be_center(i)) - dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0)) - enddo - - if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then - - continue - - else - - NAI_pol_mult_cgtos = V_n_e_cgtos(power_A(1), power_A(2), power_A(3), & - power_B(1), power_B(2), power_B(3), & - alpha, beta) - return - - endif - - p = alpha + beta - p_inv = (1.d0, 0.d0) / p - rho = alpha * beta * p_inv - - dist = (Ae_center(1) - Be_center(1)) * (Ae_center(1) - Be_center(1)) & - + (Ae_center(2) - Be_center(2)) * (Ae_center(2) - Be_center(2)) & - + (Ae_center(3) - Be_center(3)) * (Ae_center(3) - Be_center(3)) - - const_factor = dist * rho - if(real(const_factor) > 80.d0) then - NAI_pol_mult_cgtos = (0.d0, 0.d0) - return - endif - - P_center(1) = (alpha * Ae_center(1) + beta * Be_center(1)) * p_inv - P_center(2) = (alpha * Ae_center(2) + beta * Be_center(2)) * p_inv - P_center(3) = (alpha * Ae_center(3) + beta * Be_center(3)) * p_inv - - dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & - + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & - + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) - - const = p * dist_integral - factor = zexp(-const_factor) - coeff = dtwo_pi * factor * p_inv - - n_pt = 2 * ((power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3))) - if(n_pt == 0) then - NAI_pol_mult_cgtos = coeff * crint_1(0, const) - return - endif - - d(0:n_pt_in) = (0.d0, 0.d0) - call give_cpolynomial_mult_center_one_e(Ap_center, Bp_center, alpha, beta, & - power_A, power_B, C_center, n_pt_in, d, n_pt_out) - - if(n_pt_out < 0) then - NAI_pol_mult_cgtos = (0.d0, 0.d0) - return - endif - - NAI_pol_mult_cgtos = coeff * crint_sum(n_pt_out, const, d) - - return -end - -! --- - -subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & - power_A, power_B, C_center, n_pt_in, d, n_pt_out) - - BEGIN_DOC - ! Returns the explicit polynomial in terms of the "t" variable of the following - ! - ! $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, power_A(3), power_B(3) - double precision, intent(in) :: C_center(3) - complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3) - integer, intent(out) :: n_pt_out - complex*16, intent(inout) :: 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, n_pt_tmp - complex*16 :: p, P_center(3), rho, p_inv, p_inv_2 - complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2) - complex*16 :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in) - - ASSERT (n_pt_in > 1) - - p = alpha + beta - p_inv = (1.d0, 0.d0) / p - p_inv_2 = 0.5d0 * p_inv - - P_center(1) = (alpha * A_center(1) + beta * B_center(1)) * p_inv - P_center(2) = (alpha * A_center(2) + beta * B_center(2)) * p_inv - P_center(3) = (alpha * A_center(3) + beta * B_center(3)) * p_inv - - do i = 0, n_pt_in - d(i) = (0.d0, 0.d0) - d1(i) = (0.d0, 0.d0) - d2(i) = (0.d0, 0.d0) - d3(i) = (0.d0, 0.d0) - enddo - - ! --- - - n_pt1 = n_pt_in - - R1x(0) = (P_center(1) - A_center(1)) - R1x(1) = (0.d0, 0.d0) - R1x(2) = -(P_center(1) - C_center(1)) - - R1xp(0) = (P_center(1) - B_center(1)) - R1xp(1) = (0.d0, 0.d0) - R1xp(2) = -(P_center(1) - C_center(1)) - - R2x(0) = p_inv_2 - R2x(1) = (0.d0, 0.d0) - R2x(2) = -p_inv_2 - - a_x = power_A(1) - b_x = power_B(1) - - call I_x1_pol_mult_one_e_cgtos(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 - d(i) = (0.d0, 0.d0) - enddo - return - endif - - ! --- - - n_pt2 = n_pt_in - - R1x(0) = (P_center(2) - A_center(2)) - R1x(1) = (0.d0, 0.d0) - R1x(2) = -(P_center(2) - C_center(2)) - - R1xp(0) = (P_center(2) - B_center(2)) - R1xp(1) = (0.d0, 0.d0) - R1xp(2) = -(P_center(2) - C_center(2)) - - a_y = power_A(2) - b_y = power_B(2) - - call I_x1_pol_mult_one_e_cgtos(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 - d(i) = (0.d0, 0.d0) - enddo - return - endif - - ! --- - - n_pt3 = n_pt_in - - R1x(0) = (P_center(3) - A_center(3)) - R1x(1) = (0.d0, 0.d0) - R1x(2) = -(P_center(3) - C_center(3)) - - R1xp(0) = (P_center(3) - B_center(3)) - R1xp(1) = (0.d0, 0.d0) - R1xp(2) = -(P_center(3) - C_center(3)) - - a_z = power_A(3) - b_z = power_B(3) - - call I_x1_pol_mult_one_e_cgtos(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, 0.d0) - enddo - return - endif - - ! --- - - n_pt_tmp = 0 - call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp) - - n_pt_out = 0 - d1(0:n_pt_tmp) = (0.d0, 0.d0) - call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out) - do i = 0, n_pt_out - d(i) = d1(i) - enddo - -end - -! --- - -recursive subroutine I_x1_pol_mult_one_e_cgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in) - - BEGIN_DOC - ! Recursive routine involved in the electron-nucleus potential - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: a, c, n_pt_in - complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:n_pt_in) - - integer :: nx, ix, dim, iy, ny - complex*16 :: X(0:max_dim) - complex*16 :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - - dim = n_pt_in - - if((a == 0) .and. (c == 0)) then - - nd = 0 - d(0) = (1.d0, 0.d0) - return - - elseif((c < 0) .or. (nd < 0)) then - - nd = -1 - return - - elseif((a == 0) .and. (c .ne. 0)) then - - call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in) - - elseif(a == 1) then - - nx = nd - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - Y(ix) = (0.d0, 0.d0) - enddo - - call I_x2_pol_mult_one_e_cgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in) - - do ix = 0, nx - X(ix) *= dble(c) - enddo - - call multiply_cpoly(X, nx, R2x, 2, d, nd) - - ny = 0 - call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in) - call multiply_cpoly(Y, ny, R1x, 2, d, nd) - - else - - nx = 0 - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - Y(ix) = (0.d0, 0.d0) - enddo - - call I_x1_pol_mult_one_e_cgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in) - - do ix = 0, nx - X(ix) *= dble(a-1) - enddo - call multiply_cpoly(X, nx, R2x, 2, d, nd) - - nx = nd - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - - call I_x1_pol_mult_one_e_cgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in) - do ix = 0, nx - X(ix) *= dble(c) - enddo - - call multiply_cpoly(X, nx, R2x, 2, d, nd) - - ny = 0 - call I_x1_pol_mult_one_e_cgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in) - call multiply_cpoly(Y, ny, R1x, 2, d, nd) - - endif - -end - -! --- - -recursive subroutine I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, dim) - - BEGIN_DOC - ! Recursive routine involved in the electron-nucleus potential - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: dim, c - complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2) - integer, intent(inout) :: nd - complex*16, intent(out) :: d(0:max_dim) - - integer :: i, nx, ix, ny - complex*16 :: X(0:max_dim), Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - - if(c == 0) then - - nd = 0 - d(0) = (1.d0, 0.d0) - return - - elseif((nd < 0) .or. (c < 0)) then - - nd = -1 - return - - else - - nx = 0 - do ix = 0, dim - X(ix) = (0.d0, 0.d0) - Y(ix) = (0.d0, 0.d0) - enddo - - call I_x1_pol_mult_one_e_cgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim) - - do ix = 0, nx - X(ix) *= dble(c-1) - enddo - - call multiply_cpoly(X, nx, R2x, 2, d, nd) - - ny = 0 - do ix = 0, dim - Y(ix) = (0.d0, 0.d0) - enddo - - call I_x1_pol_mult_one_e_cgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim) - - if(ny .ge. 0) then - call multiply_cpoly(Y, ny, R1xp, 2, d, nd) - endif - - endif - -end - -! --- - -complex*16 function V_n_e_cgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta) - - BEGIN_DOC - ! Primitve nuclear attraction between the two primitves centered on the same atom. - ! - ! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$ - ! - ! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$ - END_DOC - - implicit none - - integer, intent(in) :: a_x, a_y, a_z, b_x, b_y, b_z - complex*16, intent(in) :: alpha, beta - - double precision :: V_phi, V_theta - complex*16 :: V_r_cgtos - - if( (iand(a_x + b_x, 1) == 1) .or. & - (iand(a_y + b_y, 1) == 1) .or. & - (iand(a_z + b_z, 1) == 1) ) then - - V_n_e_cgtos = (0.d0, 0.d0) - - else - - V_n_e_cgtos = V_r_cgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) & - * V_phi(a_x + b_x, a_y + b_y) & - * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1) - endif - -end - -! --- - -complex*16 function V_r_cgtos(n, alpha) - - BEGIN_DOC - ! Computes the radial part of the nuclear attraction integral: - ! - ! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$ - ! - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer , intent(in) :: n - complex*16, intent(in) :: alpha - - double precision :: fact - - if(iand(n, 1) .eq. 1) then - V_r_cgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1)) - else - V_r_cgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1) - endif - -end - -! --- - diff --git a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f deleted file mode 100644 index f2557fc3..00000000 --- a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f +++ /dev/null @@ -1,318 +0,0 @@ - -! --- - - BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_x, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_y, (ao_num, ao_num)] -&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_z, (ao_num, ao_num)] - - implicit none - integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) - double precision :: c, deriv_tmp - double precision :: KA2, phiA - double precision :: KB2, phiB - double precision :: aa - complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1 - complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2 - complex*16 :: xa - complex*16 :: overlap_x, overlap_y, overlap_z, overlap - complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1 - complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2 - complex*16 :: overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1 - complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2 - complex*16 :: deriv_tmp_1, deriv_tmp_2 - - - dim1 = 100 - - ! -- Dummy call to provide everything - - Ae_center(:) = (0.0d0, 0.d0) - Be_center(:) = (1.0d0, 0.d0) - Ap_center(:) = (0.0d0, 0.d0) - Bp_center(:) = (1.0d0, 0.d0) - alpha = (1.0d0, 0.d0) - beta = (0.1d0, 0.d0) - power_A = 1 - power_B = 0 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) - - ! --- - - !$OMP PARALLEL & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, & - !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, & - !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, & - !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & - !$OMP overlap_x, overlap_y, overlap_z, overlap, & - !$OMP overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1, & - !$OMP overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2, & - !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & - !$OMP overlap_y0_2, overlap_z0_2) & - !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & - !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) - !$OMP DO SCHEDULE(GUIDED) - do j = 1, ao_num - - jj = ao_nucl(j) - power_A(1) = ao_power(j,1) - power_A(2) = ao_power(j,2) - power_A(3) = ao_power(j,3) - - do i = 1, ao_num - - ii = ao_nucl(i) - power_B(1) = ao_power(i,1) - power_B(2) = ao_power(i,2) - power_B(3) = ao_power(i,3) - - ao_deriv2_cgtos_x(i,j) = 0.d0 - ao_deriv2_cgtos_y(i,j) = 0.d0 - ao_deriv2_cgtos_z(i,j) = 0.d0 - - do n = 1, ao_prim_num(j) - - alpha = ao_expo_cgtos_ord_transp(n,j) - alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 - Ap_center(m) = nucl_coord(jj,m) - Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) - enddo - phiA = ao_expo_phase_ord_transp(4,n,j) - KA2 = ao_expo_pw_ord_transp(4,n,j) - - do l = 1, ao_prim_num(i) - - beta = ao_expo_cgtos_ord_transp(l,i) - beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 - Bp_center(m) = nucl_coord(ii,m) - Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) - enddo - phiB = ao_expo_phase_ord_transp(4,l,i) - KB2 = ao_expo_pw_ord_transp(4,l,i) - - c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i) - - C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) - C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2)) - - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) - - ! --- - - power_A(1) = power_A(1) - 1 - if(power_A(1) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_m1_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1) - power_A(1) = power_A(1) - 1 - if(power_A(1) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_m2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_m2_2, overlap_y, overlap_z, overlap, dim1) - else - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(1) = power_A(1) + 1 - else - overlap_m1_1 = (0.d0, 0.d0) - overlap_m1_2 = (0.d0, 0.d0) - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(1) = power_A(1) + 1 - - power_A(1) = power_A(1) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_p1_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1) - power_A(1) = power_A(1) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_p2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_p2_2, overlap_y, overlap_z, overlap, dim1) - power_A(1) = power_A(1) - 2 - - aa = dble(power_A(1)) - xa = Ap_center(1) - Ae_center(1) - - deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) - deriv_tmp_1 = deriv_tmp_1 * overlap_y0_1 * overlap_z0_1 - - deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_2 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) - deriv_tmp_2 = deriv_tmp_2 * overlap_y0_2 * overlap_z0_2 - - deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) - - ao_deriv2_cgtos_x(i,j) += c * deriv_tmp - - ! --- - - power_A(2) = power_A(2) - 1 - if(power_A(2) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_m1_1, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1) - power_A(2) = power_A(2) - 1 - if(power_A(2) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_m2_1, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_m2_2, overlap_z, overlap, dim1) - else - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(2) = power_A(2) + 1 - else - overlap_m1_1 = (0.d0, 0.d0) - overlap_m1_2 = (0.d0, 0.d0) - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(2) = power_A(2) + 1 - - power_A(2) = power_A(2) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_p1_1, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1) - power_A(2) = power_A(2) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_p2_1, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_p2_2, overlap_z, overlap, dim1) - power_A(2) = power_A(2) - 2 - - aa = dble(power_A(2)) - xa = Ap_center(2) - Ae_center(2) - - deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) - deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_z0_1 - - deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_2 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) - deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_z0_2 - - deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) - - ao_deriv2_cgtos_y(i,j) += c * deriv_tmp - - ! --- - - power_A(3) = power_A(3) - 1 - if(power_A(3) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_y, overlap_m1_1, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1) - power_A(3) = power_A(3) - 1 - if(power_A(3) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_y, overlap_m2_1, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m2_2, overlap, dim1) - else - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(3) = power_A(3) + 1 - else - overlap_m1_1 = (0.d0, 0.d0) - overlap_m1_2 = (0.d0, 0.d0) - overlap_m2_1 = (0.d0, 0.d0) - overlap_m2_2 = (0.d0, 0.d0) - endif - power_A(3) = power_A(3) + 1 - - power_A(3) = power_A(3) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_y, overlap_p1_1, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1) - power_A(3) = power_A(3) + 1 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & - Ap_center, Bp_center, overlap_x, overlap_y, overlap_p2_1, overlap, dim1) - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & - Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p2_2, overlap, dim1) - power_A(3) = power_A(3) - 2 - - aa = dble(power_A(3)) - xa = Ap_center(3) - Ae_center(3) - - deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) - deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_y0_1 - - deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & - + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_2 & - + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) - deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_y0_2 - - deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) - - ao_deriv2_cgtos_z(i,j) += c * deriv_tmp - - ! --- - - enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cgtos, (ao_num, ao_num)] - - BEGIN_DOC - ! - ! Kinetic energy integrals in the cgtos |AO| basis. - ! - ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ - ! - END_DOC - - implicit none - - integer :: i, j - - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP PRIVATE(i, j) & - !$OMP SHARED(ao_num, ao_kinetic_integrals_cgtos, ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) - do j = 1, ao_num - do i = 1, ao_num - ao_kinetic_integrals_cgtos(i,j) = -0.5d0 * (ao_deriv2_cgtos_x(i,j) + & - ao_deriv2_cgtos_y(i,j) + & - ao_deriv2_cgtos_z(i,j)) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - -! --- - 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 dddf98d4..bd25cd76 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 @@ -9,13 +9,11 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) END_DOC double precision, intent(in) :: mu_in,C_center(3) double precision, intent(out) :: integrals_ao(ao_num,ao_num) - double precision :: NAI_pol_mult_erf_ao integer :: i,j,l,k,m - do k = 1, ao_num - do m = 1, ao_num - integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center) - enddo - enddo + double precision, allocatable :: integrals_ao_cart(:,:) + allocate(integrals_ao_cart(ao_cart_num,ao_cart_num)) + call give_all_erf_kl_ao_cart(integrals_ao_cart,mu_in,C_center) + call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num) end ! --- @@ -33,35 +31,39 @@ double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) 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 - - double precision :: NAI_pol_mult_erf - - 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 - - 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) - enddo - enddo - + double precision, allocatable :: integrals_ao(:,:),integrals_ao_cart(:,:) + allocate(integrals_ao(ao_num,ao_num),integrals_ao_cart(ao_cart_num,ao_cart_num)) + call give_all_erf_kl_ao_cart(integrals_ao_cart,mu_in,C_center) + call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num) + NAI_pol_mult_erf_ao = integrals_ao(i_ao, j_ao) end function NAI_pol_mult_erf_ao ! --- +subroutine all_NAI_pol_mult_erf_ao_with1s(beta, B_center, mu_in, C_center) + + BEGIN_DOC + ! + ! Computes ALL 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 + double precision, intent(in) :: beta, B_center(3) + double precision, intent(in) :: mu_in, C_center(3) + double precision, intent(out) :: integrals_ao(ao_num,ao_num) + + double precision :: NAI_pol_mult_erf_ao_cart_with1s + integer :: i,j + double precision, allocatable :: integrals_ao_cart(:,:) + allocate(integrals_ao_cart(ao_cart_num,ao_cart_num)) + do i = 1, ao_cart_num + do j = 1, ao_cart_num + integrals_ao_car(j,i) = NAI_pol_mult_erf_ao_cart_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center) + enddo + enddo + call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num) +end double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center) @@ -76,719 +78,17 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, 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, coef1, integral - - double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao - - ASSERT(beta .ge. 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) - coef1 = ao_coef_normalized_ordered_transp(i,i_ao) - - do j = 1, ao_prim_num(j_ao) - alpha2 = ao_expo_ordered_transp(j,j_ao) - coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao) - if(dabs(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 + double precision :: NAI_pol_mult_erf_ao_cart_with1s + integer :: i,j + double precision, allocatable :: integrals_ao(:,:),integrals_ao_cart(:,:) + allocate(integrals_ao(ao_num,ao_num),integrals_ao_cart(ao_cart_num,ao_cart_num)) + do i = 1, ao_cart_num + do j = 1, ao_cart_num + integrals_ao_car(j,i) = NAI_pol_mult_erf_ao_cart_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center) + enddo enddo + call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num) + NAI_pol_mult_erf_ao_with1s = integrals_ao(i_ao, j_ao) 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 |}$. - ! - END_DOC - - include 'utils/constants.include.F' - - 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)) - enddo - const_factor = dist * rho - if(const_factor > 80.d0) then - NAI_pol_mult_erf = 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_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 - ! 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, 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 - NAI_pol_mult_erf = 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 = accu * coeff - -end function NAI_pol_mult_erf - -! --- - -subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points) - - 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 |}$. - ! - END_DOC - - include 'utils/constants.include.F' - - implicit none - - integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv - integer, intent(in) :: power_A(3), power_B(3) - double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in - double precision, intent(in) :: C_center(LD_C,3) - double precision, intent(out) :: res_v(LD_resv) - - integer :: i, n_pt, n_pt_out, ipoint - 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, p_new2, coef_tmp - - double precision :: rint - - res_V(1:LD_resv) = 0.d0 - - p = alpha + beta - p_inv = 1.d0 / p - p_inv_2 = 0.5d0 * p_inv - rho = alpha * beta * p_inv - p_new = mu_in / dsqrt(p + mu_in * mu_in) - p_new2 = p_new * p_new - coef_tmp = p * p_new2 - - dist = 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)) - enddo - - const_factor = dist * rho - if(const_factor > 80.d0) then - return - endif - 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) ) - - if(n_pt == 0) then - - do ipoint = 1, n_points - dist_integral = 0.d0 - do i = 1, 3 - dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i)) - enddo - const = coef_tmp * dist_integral - - res_v(ipoint) = coeff * rint(0, const) - enddo - - else - - do ipoint = 1, n_points - dist_integral = 0.d0 - do i = 1, 3 - dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i)) - enddo - const = coef_tmp * dist_integral - - do i = 0, n_pt_in - d(i) = 0.d0 - enddo - call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center) - - if(n_pt_out < 0) then - res_v(ipoint) = 0.d0 - cycle - 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 - - res_v(ipoint) = accu * coeff - enddo - - endif - -end subroutine NAI_pol_mult_erf_v - -! --- - -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 - - - ! 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 - A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv - A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv - A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv - dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) & - + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) & - + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) - - 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 - P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv - P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv - P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv - dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) & - + (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) & - + (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) - - const_factor = const_factor12 + dist * rho - if(const_factor > 80.d0) then - NAI_pol_mult_erf_with1s = 0.d0 - return - endif - - dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & - + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & - + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) - - ! --- - - 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 NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points) - - 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, LD_B, LD_C, LD_resv, n_points - integer, intent(in) :: power_A1(3), power_A2(3) - double precision, intent(in) :: A1_center(3), A2_center(3) - double precision, intent(in) :: C_center(LD_C,3), B_center(LD_B,3) - double precision, intent(in) :: alpha1, alpha2, beta, mu_in - double precision, intent(out) :: res_v(LD_resv) - - integer :: i, n_pt, n_pt_out, ipoint - 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, p_new2, coef_tmp, cons_tmp - - double precision :: rint - - - res_V(1:LD_resv) = 0.d0 - - ! 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 - A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv - A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv - A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv - dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1))& - + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2))& - + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) - - const_factor12 = dist12 * rho12 - if(const_factor12 > 80.d0) then - 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 - p_new = mu_in / dsqrt(p + mu_in * mu_in) - p_new2 = p_new * p_new - coef_tmp = dtwo_pi * p_inv * p_new - cons_tmp = p * p_new2 - n_pt = 2 * (power_A1(1) + power_A2(1) + power_A1(2) + power_A2(2) + power_A1(3) + power_A2(3) ) - - if(n_pt == 0) then - - do ipoint = 1, n_points - - dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))& - + (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))& - + (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3)) - const_factor = const_factor12 + dist * rho - if(const_factor > 80.d0) cycle - coeff = coef_tmp * dexp(-const_factor) - - P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv - P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv - P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv - dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))& - + (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))& - + (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3)) - const = cons_tmp * dist_integral - - res_v(ipoint) = coeff * rint(0, const) - enddo - - else - - do ipoint = 1, n_points - - dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))& - + (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))& - + (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3)) - const_factor = const_factor12 + dist * rho - if(const_factor > 80.d0) cycle - coeff = coef_tmp * dexp(-const_factor) - - P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv - P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv - P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv - dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))& - + (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))& - + (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3)) - const = cons_tmp * dist_integral - - do i = 0, n_pt_in - d(i) = 0.d0 - enddo - !TODO: VECTORIZE HERE - call give_polynomial_mult_center_one_e_erf_opt(A1_center, A2_center, power_A1, power_A2, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center) - - if(n_pt_out < 0) then - cycle - 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 - - res_v(ipoint) = accu * coeff - enddo - - endif - -end subroutine NAI_pol_mult_erf_with1s_v - -! --- - -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(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) - - 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 = (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.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 - - do i = 0, n_pt_in - d (i) = 0.d0 - d1(i) = 0.d0 - d2(i) = 0.d0 - d3(i) = 0.d0 - enddo - - 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) - if(n_pt1<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - 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 = (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 = (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) - if(n_pt2<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - 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 = (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 - !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 - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - - n_pt_tmp = 0 - 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) - do i = 0, n_pt_out - d(i) = d1(i) - enddo - -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) - - 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) - 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, pq_inv, p10_1, p10_2, p01_1, p01_2 - double precision :: p,P_center(3),rho,p_inv,p_inv_2 - accu = 0.d0 - ASSERT (n_pt_in > 1) - p = alpha+beta - p_inv = 1.d0/p - p_inv_2 = 0.5d0/p - do i =1, 3 - P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv - enddo - - 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))* mu_in**2 / (p+mu_in*mu_in) - ! 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))* mu_in**2 / (p+mu_in*mu_in) - !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* mu_in**2 / (p+mu_in*mu_in) - !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 - 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) - ! print*,'passed the first I_x1' - if(n_pt1<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - - R1x(0) = (P_center(2) - A_center(2)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in) - ! 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))* mu_in**2 / (p+mu_in*mu_in) - !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) - ! print*,'passed the second I_x1' - if(n_pt2<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - - - R1x(0) = (P_center(3) - A_center(3)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in) - ! 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))* mu_in**2 / (p+mu_in*mu_in) - !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 - a_z = power_A(3) - b_z = power_B(3) - - ! print*,'a_z = ',a_z - ! print*,'b_z = ',b_z - call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) - ! print*,'passed the third I_x1' - 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 - d1(i) = 0.d0 - enddo - n_pt_out = 0 - 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 - diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index b4b4aa02..df04cac9 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -27,78 +27,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] else - if(use_cgtos) then - - do j = 1, ao_num - do i = 1, ao_num - ao_integrals_n_e(i,j) = ao_integrals_n_e_cgtos(i,j) - enddo - enddo - - else - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c,c1,n_pt_in) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& - !$OMP n_pt_max_integrals,ao_integrals_n_e,nucl_num,nucl_charge) - - n_pt_in = n_pt_max_integrals - - !$OMP DO SCHEDULE (dynamic) - - do j = 1, ao_num - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - 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) - - double precision :: c, c1 - c = 0.d0 - - do k = 1, nucl_num - double precision :: Z - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - !print *, ' ' - !print *, A_center, B_center, C_center, power_A, power_B - !print *, alpha, beta - - c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & - , alpha, beta, C_center, n_pt_in ) - - !print *, ' c1 = ', c1 - - c = c - Z * c1 - - enddo - ao_integrals_n_e(i,j) = ao_integrals_n_e(i,j) & - + ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) * c - enddo - enddo - enddo - enddo - - !$OMP END DO - !$OMP END PARALLEL - - endif - + call ao_cart_to_ao_basis(ao_cart_coul_n_e, ao_cart_num, ao_integrals_n_e, ao_num) IF(do_pseudo) THEN ao_integrals_n_e += ao_pseudo_integrals @@ -147,487 +76,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nuc ! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle` END_DOC implicit none - double precision :: alpha, beta - integer :: i_c,num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m - double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - - ao_integrals_n_e_per_atom = 0.d0 - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,power_A,power_B,& - !$OMP num_A,num_B,c,n_pt_in,C_center) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& - !$OMP n_pt_max_integrals,ao_integrals_n_e_per_atom,nucl_num) - n_pt_in = n_pt_max_integrals - !$OMP DO SCHEDULE (dynamic) - - double precision :: c - do j = 1, ao_num - power_A(1)= ao_power(j,1) - power_A(2)= ao_power(j,2) - power_A(3)= ao_power(j,3) - num_A = ao_nucl(j) - A_center(1) = nucl_coord(num_A,1) - A_center(2) = nucl_coord(num_A,2) - A_center(3) = nucl_coord(num_A,3) - do k = 1, nucl_num - C_center(1) = nucl_coord(k,1) - C_center(2) = nucl_coord(k,2) - C_center(3) = nucl_coord(k,3) - do i = 1, ao_num - power_B(1)= ao_power(i,1) - power_B(2)= ao_power(i,2) - power_B(3)= ao_power(i,3) - num_B = ao_nucl(i) - B_center(1) = nucl_coord(num_B,1) - B_center(2) = nucl_coord(num_B,2) - B_center(3) = nucl_coord(num_B,3) - c = 0.d0 - 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) - c = c + NAI_pol_mult(A_center,B_center,power_A,power_B, & - alpha,beta,C_center,n_pt_in) & - * ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) - enddo - enddo - ao_integrals_n_e_per_atom(i,j,k) = -c - enddo - enddo + integer :: i + do i = 1, nucl_num + call ao_cart_to_ao_basis(ao_cart_coul_n_e_per_atom(1,1,i), ao_cart_num,ao_integrals_n_e_per_atom(1,1,i), ao_num) enddo - !$OMP END DO - !$OMP END PARALLEL - END_PROVIDER - - -double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - BEGIN_DOC -! Computes the electron-nucleus attraction with two primitves. -! -! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle` - 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 - integer :: 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,rho,dist,const,pouet_2,p,p_inv,factor - double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi - double precision :: V_n_e,const_factor,dist_integral,tmp - double precision :: accu,epsilo,rint - integer :: n_pt_out,lmax - include 'utils/constants.include.F' - if ( (A_center(1)/=B_center(1)).or. & - (A_center(2)/=B_center(2)).or. & - (A_center(3)/=B_center(3)).or. & - (A_center(1)/=C_center(1)).or. & - (A_center(2)/=C_center(2)).or. & - (A_center(3)/=C_center(3))) then - continue - else - NAI_pol_mult = V_n_e(power_A(1),power_A(2),power_A(3), & - power_B(1),power_B(2),power_B(3),alpha,beta) - return - endif - p = alpha + beta - p_inv = 1.d0/p - 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)) - enddo - const_factor = dist*rho - const = p * dist_integral - if(const_factor > 80.d0)then - NAI_pol_mult = 0.d0 - return - endif - factor = dexp(-const_factor) - coeff = dtwo_pi * factor * p_inv - lmax = 20 - - ! print*, "b" - 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)) ) - if (n_pt == 0) then - epsilo = 1.d0 - pouet = rint(0,const) - NAI_pol_mult = coeff * pouet - return - endif - - call give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out) - - - if(n_pt_out<0)then - NAI_pol_mult = 0.d0 - return - endif - accu = 0.d0 - - ! 1/r1 standard attraction integral - epsilo = 1.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) - enddo - NAI_pol_mult = accu * coeff - -end - -! --- - -subroutine give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out) - implicit none - BEGIN_DOC - ! Returns the explicit polynomial in terms of the "t" variable of the following - ! - ! $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 - 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) - double precision, intent(in) :: alpha,beta - 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, pq_inv, p10_1, p10_2, p01_1, p01_2 - double precision :: p,P_center(3),rho,p_inv,p_inv_2 - - accu = 0.d0 - - ASSERT (n_pt_in > 1) - p = alpha+beta - p_inv = 1.d0/p - p_inv_2 = 0.5d0/p - do i =1, 3 - P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv - enddo - - 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)) - - R1xp(0) = (P_center(1) - B_center(1)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(1) - C_center(1)) - - R2x(0) = p_inv_2 - R2x(1) = 0.d0 - R2x(2) = -p_inv_2 - - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - do i = 0,n_pt_in - 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) - - if(n_pt1<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - - R1x(0) = (P_center(2) - A_center(2)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(2) - C_center(2)) - - R1xp(0) = (P_center(2) - B_center(2)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(2) - C_center(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) - - if(n_pt2<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - - - R1x(0) = (P_center(3) - A_center(3)) - R1x(1) = 0.d0 - R1x(2) = -(P_center(3) - C_center(3)) - - R1xp(0) = (P_center(3) - B_center(3)) - R1xp(1) = 0.d0 - R1xp(2) =-(P_center(3) - C_center(3)) - - 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 - 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 - d1(i) = 0.d0 - enddo - n_pt_out = 0 - 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 - - -recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) - implicit none - BEGIN_DOC -! Recursive routine involved in the electron-nucleus potential - END_DOC - integer , intent(in) :: n_pt_in - double precision,intent(inout) :: d(0:n_pt_in) - integer,intent(inout) :: nd - integer, intent(in) :: a,c - double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2) - include 'utils/constants.include.F' - double precision :: X(0:max_dim) - double precision :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - integer :: nx, ix,dim,iy,ny - dim = n_pt_in - ! print*,'a,c = ',a,c - ! print*,'nd_in = ',nd - - if( (a==0) .and. (c==0))then - nd = 0 - d(0) = 1.d0 - return - elseif( (c<0).or.(nd<0) )then - nd = -1 - return - else if ((a==0).and.(c.ne.0)) then - call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,n_pt_in) - else if (a==1) then - nx = nd - do ix=0,n_pt_in - X(ix) = 0.d0 - Y(ix) = 0.d0 - enddo - call I_x2_pol_mult_one_e(c-1,R1x,R1xp,R2x,X,nx,n_pt_in) - do ix=0,nx - X(ix) *= dble(c) - enddo -! call multiply_poly(X,nx,R2x,2,d,nd) - call multiply_poly_c2(X,nx,R2x,d,nd) - ny=0 - call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in) -! call multiply_poly(Y,ny,R1x,2,d,nd) - call multiply_poly_c2(Y,ny,R1x,d,nd) - else - do ix=0,n_pt_in - X(ix) = 0.d0 - Y(ix) = 0.d0 - enddo - nx = 0 - call I_x1_pol_mult_one_e(a-2,c,R1x,R1xp,R2x,X,nx,n_pt_in) - do ix=0,nx - X(ix) *= dble(a-1) - enddo -! call multiply_poly(X,nx,R2x,2,d,nd) - call multiply_poly_c2(X,nx,R2x,d,nd) - - nx = nd - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - call I_x1_pol_mult_one_e(a-1,c-1,R1x,R1xp,R2x,X,nx,n_pt_in) - do ix=0,nx - X(ix) *= dble(c) - enddo -! call multiply_poly(X,nx,R2x,2,d,nd) - call multiply_poly_c2(X,nx,R2x,d,nd) - ny=0 - call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in) -! call multiply_poly(Y,ny,R1x,2,d,nd) - call multiply_poly_c2(Y,ny,R1x,d,nd) - endif -end - -recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim) - implicit none - BEGIN_DOC -! Recursive routine involved in the electron-nucleus potential - END_DOC - integer , intent(in) :: dim - include 'utils/constants.include.F' - double precision :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: c - double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2) - integer :: i - - if(c==0) then - nd = 0 - d(0) = 1.d0 - return - elseif ((nd<0).or.(c<0))then - nd = -1 - return - else - integer :: nx, ix,ny - double precision :: X(0:max_dim),Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - do ix=0,dim - X(ix) = 0.d0 - Y(ix) = 0.d0 - enddo - nx = 0 - call I_x1_pol_mult_one_e(0,c-2,R1x,R1xp,R2x,X,nx,dim) - do ix=0,nx - X(ix) *= dble(c-1) - enddo -! call multiply_poly(X,nx,R2x,2,d,nd) - call multiply_poly_c2(X,nx,R2x,d,nd) - ny = 0 - do ix=0,dim - Y(ix) = 0.d0 - enddo - - call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim) - if(ny.ge.0)then -! call multiply_poly(Y,ny,R1xp,2,d,nd) - call multiply_poly_c2(Y,ny,R1xp,d,nd) - endif - endif -end - -double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta) - implicit none - BEGIN_DOC -! Primitve nuclear attraction between the two primitves centered on the same atom. -! -! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$ -! -! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$ - END_DOC - integer :: a_x,a_y,a_z,b_x,b_y,b_z - double precision :: alpha,beta - double precision :: V_r, V_phi, V_theta - if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then - V_n_e = 0.d0 - else - V_n_e = V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta) & - * V_phi(a_x+b_x,a_y+b_y) & - * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1) - endif - -end - - -double precision function int_gaus_pol(alpha,n) - implicit none - BEGIN_DOC -! Computes the integral: -! -! $\int_{-\infty}^{\infty} x^n \exp(-\alpha x^2) dx$. - END_DOC - double precision :: alpha - integer :: n - double precision :: dble_fact - include 'utils/constants.include.F' - - int_gaus_pol = 0.d0 - if(iand(n,1).eq.0)then - int_gaus_pol = dsqrt(alpha/pi) - double precision :: two_alpha - two_alpha = alpha+alpha - integer :: i - do i=1,n,2 - int_gaus_pol = int_gaus_pol * two_alpha - enddo - int_gaus_pol = dble_fact(n -1) / int_gaus_pol - endif - -end - -double precision function V_r(n,alpha) - implicit none - BEGIN_DOC - ! Computes the radial part of the nuclear attraction integral: - ! - ! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$ - ! - END_DOC - double precision :: alpha, fact - integer :: n - include 'utils/constants.include.F' - if(iand(n,1).eq.1)then - V_r = 0.5d0 * fact(shiftr(n,1)) / (alpha ** (shiftr(n,1) + 1)) - else - V_r = sqpi * fact(n) / fact(shiftr(n,1)) * (0.5d0/sqrt(alpha)) ** (n+1) - endif -end - - - -BEGIN_PROVIDER [ double precision, ao_sphe_integrals_n_e, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! |AO| VneVne inntegrals matrix in the spherical basis set - END_DOC - double precision, allocatable :: tmp(:,:) - allocate (tmp(ao_sphe_num,ao_num)) - - call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & - ao_integrals_n_e,size(ao_integrals_n_e,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - tmp, size(tmp,1), & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & - ao_sphe_integrals_n_e,size(ao_sphe_integrals_n_e,1)) - - deallocate(tmp) - -END_PROVIDER diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index f752614c..77852359 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -8,16 +8,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)] call ezfio_get_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals) print *, 'AO pseudopotential integrals read from disk' else - - ao_pseudo_integrals = 0.d0 - if (do_pseudo) then - if (pseudo_klocmax > 0) then - ao_pseudo_integrals += ao_pseudo_integrals_local - endif - if (pseudo_kmax > 0) then - ao_pseudo_integrals += ao_pseudo_integrals_non_local - endif - endif + call ao_cart_to_ao_basis(ao_cart_pseudo_integrals, ao_cart_num, ao_pseudo_integrals, ao_num) endif if (write_ao_integrals_pseudo) then @@ -27,349 +18,3 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)] END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] - use omp_lib - implicit none - BEGIN_DOC - ! Local pseudo-potential - END_DOC - include 'utils/constants.include.F' - double precision :: alpha, beta, gama, delta - integer :: num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,m - double precision :: Vloc, Vpseudo - - double precision :: wall_1, wall_2, wall_0 - integer :: thread_num - double precision :: c - double precision :: Z - - PROVIDE ao_coef_normalized_ordered_transp - PROVIDE pseudo_v_k_transp pseudo_n_k_transp pseudo_klocmax pseudo_dz_k_transp - - ao_pseudo_integrals_local = 0.d0 - - print*, 'Providing the nuclear electron pseudo integrals (local)' - - ! Dummy iteration for OpenMP - j=1 - i=1 - l=1 - m=1 - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - alpha = ao_expo_ordered_transp(l,j) - beta = ao_expo_ordered_transp(m,i) - c = 0.d0 - do k = 1, nucl_num - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - c = c + Vloc(pseudo_klocmax, & - pseudo_v_k_transp (1,k), & - pseudo_n_k_transp (1,k), & - pseudo_dz_k_transp(1,k), & - A_center,power_A,alpha,B_center,power_B,beta,C_center) - - enddo - - - ao_pseudo_integrals_local = 0.d0 - call wall_time(wall_1) - - thread_num = 0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c, & - !$OMP wall_0,wall_2,thread_num) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& - !$OMP ao_pseudo_integrals_local,nucl_num,nucl_charge, & - !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k_transp,pseudo_n_k_transp, pseudo_dz_k_transp,& - !$OMP wall_1) - - !$ thread_num = omp_get_thread_num() - - wall_0 = wall_1 - !$OMP DO - - do j = 1, ao_num - - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - 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) - if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& - < thresh) then - cycle - endif - - beta = ao_expo_ordered_transp(m,i) - c = 0.d0 - - do k = 1, nucl_num - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - c = c + Vloc(pseudo_klocmax, & - pseudo_v_k_transp (1,k), & - pseudo_n_k_transp (1,k), & - pseudo_dz_k_transp(1,k), & - A_center,power_A,alpha,B_center,power_B,beta,C_center) - - enddo - ao_pseudo_integrals_local(i,j) = ao_pseudo_integrals_local(i,j) +& - ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c - enddo - enddo - enddo - - call wall_time(wall_2) - if (thread_num == 0) then - if (wall_2 - wall_0 > 1.d0) then - wall_0 = wall_2 - print*, 100.*float(j)/float(ao_num), '% in ', & - wall_2-wall_1, 's' - endif - endif - enddo - - !$OMP END DO - !$OMP END PARALLEL - - do i=1,ao_num - do j=1,i - ao_pseudo_integrals_local(j,i) = 0.5d0*(ao_pseudo_integrals_local(i,j) + ao_pseudo_integrals_local(i,j)) - ao_pseudo_integrals_local(i,j) = ao_pseudo_integrals_local(i,j) - enddo - enddo - END_PROVIDER - - - BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_non_local, (ao_num,ao_num)] - use omp_lib - implicit none - BEGIN_DOC - ! Non-local pseudo-potential - END_DOC - include 'utils/constants.include.F' - double precision :: alpha, beta, gama, delta - integer :: num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,m - double precision :: Vloc, Vpseudo - - double precision :: wall_1, wall_2, wall_0 - integer :: thread_num - double precision :: c - double precision :: Z - - PROVIDE ao_coef_normalized_ordered_transp - PROVIDE pseudo_lmax pseudo_kmax pseudo_v_kl_transp pseudo_n_kl_transp pseudo_dz_kl_transp - ao_pseudo_integrals_non_local = 0.d0 - - print*, 'Providing the nuclear electron pseudo integrals (non-local)' - - call wall_time(wall_1) - thread_num = 0 - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c, & - !$OMP wall_0,wall_2,thread_num) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& - !$OMP ao_pseudo_integrals_non_local,nucl_num,nucl_charge,& - !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl_transp, pseudo_v_kl_transp, pseudo_dz_kl_transp,& - !$OMP wall_1) - - !$ thread_num = omp_get_thread_num() - - wall_0 = wall_1 - !$OMP DO SCHEDULE (guided) -! - do j = 1, ao_num - - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - 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) - if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& - < thresh) then - cycle - endif - - beta = ao_expo_ordered_transp(m,i) - c = 0.d0 - - do k = 1, nucl_num - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - c = c + Vpseudo(pseudo_lmax,pseudo_kmax, & - pseudo_v_kl_transp(1,0,k), & - pseudo_n_kl_transp(1,0,k), & - pseudo_dz_kl_transp(1,0,k), & - A_center,power_A,alpha,B_center,power_B,beta,C_center) - enddo - ao_pseudo_integrals_non_local(i,j) = ao_pseudo_integrals_non_local(i,j) +& - ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c - enddo - enddo - enddo - - call wall_time(wall_2) - if (thread_num == 0) then - if (wall_2 - wall_0 > 1.d0) then - wall_0 = wall_2 - print*, 100.*float(j)/float(ao_num), '% in ', & - wall_2-wall_1, 's' - endif - endif - enddo - - !$OMP END DO - - !$OMP END PARALLEL - - - do i=1,ao_num - do j=1,i - ao_pseudo_integrals_non_local(j,i) = 0.5d0*(ao_pseudo_integrals_non_local(i,j) + ao_pseudo_integrals_non_local(i,j)) - ao_pseudo_integrals_non_local(i,j) = ao_pseudo_integrals_non_local(i,j) - enddo - enddo -END_PROVIDER - - BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ] -&BEGIN_PROVIDER [ integer , pseudo_n_k_transp, (pseudo_klocmax,nucl_num) ] -&BEGIN_PROVIDER [ double precision, pseudo_dz_k_transp, (pseudo_klocmax,nucl_num)] - implicit none - BEGIN_DOC - ! Transposed arrays for pseudopotentials - END_DOC - - integer :: i,j - do j=1,nucl_num - do i=1,pseudo_klocmax - pseudo_v_k_transp (i,j) = pseudo_v_k (j,i) - pseudo_n_k_transp (i,j) = pseudo_n_k (j,i) - pseudo_dz_k_transp(i,j) = pseudo_dz_k(j,i) - enddo - enddo -END_PROVIDER - - BEGIN_PROVIDER [ double precision, pseudo_v_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ] -&BEGIN_PROVIDER [ integer , pseudo_n_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ] -&BEGIN_PROVIDER [ double precision, pseudo_dz_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num)] - implicit none - BEGIN_DOC - ! Transposed arrays for pseudopotentials - END_DOC - - integer :: i,j,l - do j=1,nucl_num - do l=0,pseudo_lmax - do i=1,pseudo_kmax - pseudo_v_kl_transp (i,l,j) = pseudo_v_kl (j,i,l) - pseudo_n_kl_transp (i,l,j) = pseudo_n_kl (j,i,l) - pseudo_dz_kl_transp(i,l,j) = pseudo_dz_kl(j,i,l) - enddo - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_local, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! |AO| pseudo_integrals_local matrix in the spherical basis set - END_DOC - double precision, allocatable :: tmp(:,:) - allocate (tmp(ao_sphe_num,ao_num)) - - call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & - ao_pseudo_integrals_local,size(ao_pseudo_integrals_local,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - tmp, size(tmp,1), & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & - ao_sphe_pseudo_integrals_local,size(ao_sphe_pseudo_integrals_local,1)) - - deallocate(tmp) - -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_non_local, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! |AO| pseudo_integrals_non_local matrix in the spherical basis set - END_DOC - double precision, allocatable :: tmp(:,:) - allocate (tmp(ao_sphe_num,ao_num)) - - call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), & - ao_pseudo_integrals_non_local,size(ao_pseudo_integrals_non_local,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - tmp, size(tmp,1), & - ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, & - ao_sphe_pseudo_integrals_non_local,size(ao_sphe_pseudo_integrals_non_local,1)) - - deallocate(tmp) - -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals, (ao_sphe_num,ao_sphe_num)] - implicit none - BEGIN_DOC - ! Pseudo-potential integrals in the |AO| basis set. - END_DOC - - ao_sphe_pseudo_integrals = 0.d0 - if (do_pseudo) then - if (pseudo_klocmax > 0) then - ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_local - endif - if (pseudo_kmax > 0) then - ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_non_local - endif - endif - -END_PROVIDER - diff --git a/src/ao_one_e_ints/pot_pt_charges.irp.f b/src/ao_one_e_ints/pot_pt_charges.irp.f index f939b373..a0da221e 100644 --- a/src/ao_one_e_ints/pot_pt_charges.irp.f +++ b/src/ao_one_e_ints/pot_pt_charges.irp.f @@ -10,99 +10,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)] END_DOC implicit none - integer :: num_A, num_B, power_A(3), power_B(3) - integer :: i, j, k, l, n_pt_in, m - double precision :: alpha, beta - double precision :: A_center(3),B_center(3),C_center(3) - double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - - ao_integrals_pt_chrg = 0.d0 - -! if (read_ao_integrals_pt_chrg) then -! -! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) -! print *, 'AO N-e integrals read from disk' -! -! else - -! if(use_cgtos) then -! !print *, " use_cgtos for ao_integrals_pt_chrg ?", use_cgtos -! -! do j = 1, ao_num -! do i = 1, ao_num -! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cgtos(i,j) -! enddo -! enddo -! -! else - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c,c1,n_pt_in) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,& - !$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z) - - n_pt_in = n_pt_max_integrals - - !$OMP DO SCHEDULE (dynamic) - - do j = 1, ao_num - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - 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) - - double precision :: c, c1 - c = 0.d0 - - do k = 1, n_pts_charge - double precision :: Z - Z = pts_charge_z(k) - - C_center(1:3) = pts_charge_coord(k,1:3) - - c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & - , alpha, beta, C_center, n_pt_in ) - - c = c - Z * c1 - - enddo - ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) & - + ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) * c - enddo - enddo - enddo - enddo - - !$OMP END DO - !$OMP END PARALLEL - -! endif - - -! IF(do_pseudo) THEN -! ao_integrals_pt_chrg += ao_pseudo_integrals -! ENDIF - -! endif - - -! if (write_ao_integrals_pt_chrg) then -! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) -! print *, 'AO N-e integrals written to disk' -! endif + call ao_cart_to_ao_basis(ao_cart_integrals_pt_chrg, ao_cart_num, ao_integrals_pt_chrg, ao_num) END_PROVIDER diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 deleted file mode 100644 index 141292d8..00000000 --- a/src/ao_one_e_ints/pseudopot.f90 +++ /dev/null @@ -1,2127 +0,0 @@ -!! INFO : You can display equations using : http://www.codecogs.com/latex/eqneditor.php - -!! -!! {\tt Vps}(C) = \langle \Phi_A|{\tt Vloc}(C)+{\tt Vpp}(C)| \Phi_B \rangle -!! -!! with : -!! -!! {\tt Vloc}(C)=\sum_{k=1}^{\tt klocmax} v_k r_C^{n_k} \exp(-dz_k r_C^2) \\ -!! -!! {\tt Vpp}(C)=\sum_{l=0}^{\tt lmax}\left( \sum_{k=1}^{\tt kmax} v_{kl} -!! r_C^{n_{kl}} \exp(-dz_{kl} r_C)^2 \right) |l\rangle \langle l| -!! -double precision function Vps & -(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) -implicit none -integer n_a(3),n_b(3) -double precision g_a,g_b,a(3),b(3),c(3) -integer lmax,kmax,n_kl(kmax,0:lmax) -double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) -integer klocmax,n_k(klocmax) -double precision v_k(klocmax),dz_k(klocmax) -double precision Vloc,Vpseudo - -Vps=Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) & - +Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) -end -!! -!! Vps_num: brute force numerical evaluation of the same matrix element Vps -!! -double precision function Vps_num & -(npts,rmax,a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) -implicit none -integer n_a(3),n_b(3) -double precision g_a,g_b,a(3),b(3),c(3),rmax -integer lmax,kmax,n_kl(kmax,0:lmax) -double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) -integer klocmax,n_k(klocmax) -double precision v_k(klocmax),dz_k(klocmax) -double precision Vloc_num,Vpseudo_num,v1,v2 -integer npts,nptsgrid -nptsgrid=50 -call initpseudos(nptsgrid) -v1=Vloc_num(npts,rmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) -v2=Vpseudo_num(nptsgrid,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) -Vps_num=v1+v2 -end - -double precision function Vloc_num(npts_over,xmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) -implicit none -integer klocmax -double precision v_k(klocmax),dz_k(klocmax) -integer n_k(klocmax) -integer npts_over,ix,iy,iz -double precision xmax,dx,x,y,z -double precision a(3),b(3),c(3),term,r,orb_phi,g_a,g_b,ac(3),bc(3) -integer n_a(3),n_b(3),k,l -do l=1,3 - ac(l)=a(l)-c(l) - bc(l)=b(l)-c(l) -enddo -dx=2.d0*xmax/npts_over -Vloc_num=0.d0 -do ix=1,npts_over -do iy=1,npts_over -do iz=1,npts_over - x=-xmax+dx*ix+dx/2.d0 - y=-xmax+dx*iy+dx/2.d0 - z=-xmax+dx*iz+dx/2.d0 - term=orb_phi(x,y,z,n_a,ac,g_a)*orb_phi(x,y,z,n_b,bc,g_b) - r=dsqrt(x**2+y**2+z**2) - do k=1,klocmax - Vloc_num=Vloc_num+dx**3*v_k(k)*r**n_k(k)*dexp(-dz_k(k)*r**2)*term - enddo -enddo -enddo -enddo -end - -double precision function orb_phi(x,y,z,npower,center,gamma) -implicit none -integer npower(3) -double precision x,y,z,r2,gamma,center(3) -r2=(x-center(1))**2+(y-center(2))**2+(z-center(3))**2 -orb_phi=(x-center(1))**npower(1)*(y-center(2))**npower(2)*(z-center(3))**npower(3) -orb_phi=orb_phi*dexp(-gamma*r2) -end - -!! Real spherical harmonics Ylm - -! factor = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 -! Y_lm(theta,phi) = -! m > 0 factor* P_l^|m|(cos(theta)) cos (|m| phi) -! m = 0 1/sqrt(2) *factor* P_l^0(cos(theta)) -! m < 0 factor* P_l^|m|(cos(theta)) sin (|m| phi) -! -! x=cos(theta) - - double precision function ylm_real(l,m,x,phi) - implicit none - integer :: MM, iabs_m, m, l - double precision :: pi, fourpi, factor, x, phi, coef - double precision :: xchap, ychap, zchap - double precision, external :: fact - double precision :: PM(0:100,0:100), plm - MM=100 - pi=dacos(-1.d0) - fourpi=4.d0*pi - iabs_m=iabs(m) - if(iabs_m.gt.l)stop 'm must be between -l and l' - factor= dsqrt( ((l+l+1)*fact(l-iabs_m))/(fourpi*fact(l+iabs_m)) ) - if(dabs(x).gt.1.d0)then - print*,'pb. in ylm_no' - print*,'x=',x - stop - endif - call LPMN(MM,l,l,X,PM) - plm=PM(iabs_m,l) - coef=factor*plm - if(m.gt.0)ylm_real=dsqrt(2.d0)*coef*dcos(iabs_m*phi) - if(m.eq.0)ylm_real=coef - if(m.lt.0)ylm_real=dsqrt(2.d0)*coef*dsin(iabs_m*phi) - - if(l.eq.0)ylm_real=dsqrt(1.d0/fourpi) - - xchap=dsqrt(1.d0-x**2)*dcos(phi) - ychap=dsqrt(1.d0-x**2)*dsin(phi) - zchap=x - if(l.eq.1.and.m.eq.1)ylm_real=dsqrt(3.d0/fourpi)*xchap - if(l.eq.1.and.m.eq.0)ylm_real=dsqrt(3.d0/fourpi)*zchap - if(l.eq.1.and.m.eq.-1)ylm_real=dsqrt(3.d0/fourpi)*ychap - - if(l.eq.2.and.m.eq.2)ylm_real=dsqrt(15.d0/16.d0/pi)*(xchap*xchap-ychap*ychap) - if(l.eq.2.and.m.eq.1)ylm_real=dsqrt(15.d0/fourpi)*xchap*zchap - if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(2.d0*zchap*zchap-xchap*xchap-ychap*ychap) - if(l.eq.2.and.m.eq.-1)ylm_real=dsqrt(15.d0/fourpi)*ychap*zchap - if(l.eq.2.and.m.eq.-2)ylm_real=dsqrt(15.d0/fourpi)*xchap*ychap - - if(l.gt.2)stop 'l > 2 not coded!' - - end -! _ -! | | -! __ __ _ __ ___ ___ _ _ __| | ___ -! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ -! \ V / | |_) \__ \ __/ |_| | (_| | (_) | -! \_/ | .__/|___/\___|\__,_|\____|\___/ -! | | -! |_| - -!! Routine Vpseudo is based on formumla (66) -!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976): -!! -!! Vpseudo= (4pi)**2* \sum_{l=0}^lmax \sum_{m=-l}^{l} -!! \sum{lambda=0}^{l+nA} \sum_{mu=-lambda}^{lambda} -!! \sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} -!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3)* Y_{lambda mu}(AC_unit) -!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* -!! bigI(lambda,mu,l,m,k1,k2,k3) -!! \sum{lambdap=0}^{l+nB} \sum_{mup=-lambdap}^{lambdap} -!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} -!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p)* Y_{lambdap mup}(BC_unit) -!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* -!! bigI(lambdap,mup,l,m,k1p,k2p,k3p)* -!! \sum_{k=1}^{kmax} v_kl(k,l)* -!! bigR(lambda,lambdap,k1+k2+k3+k1p+k2p+k3p+n_kl(k,l),g_a,g_b,AC,BC,dz_kl(k,l)) -!! -!! nA=nAx+nAy+nAz -!! nB=nBx+nBy+nBz -!! AC=|A-C| -!! AC_x= A_x-C_x, etc. -!! BC=|B-C| -!! AC_unit= vect(AC)/AC -!! BC_unit= vect(BC)/BC -!! bigI(lambda,mu,l,m,k1,k2,k3)= -!! \int dOmega Y_{lambda mu}(Omega) xchap^k1 ychap^k2 zchap^k3 Y_{l m}(Omega) -!! -!! bigR(lambda,lambdap,N,g_a,g_b,gamm_k,AC,BC) -!! = exp(-g_a* AC**2 -g_b* BC**2) * int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) -!! /int dx x^{ktot} exp(-g_k)*x**2) M_lambda(2 g_k D x) - -double precision function Vpseudo & -(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) -implicit none - -! ___ -! | ._ ._ _|_ -! _|_ | | |_) |_| |_ -! | -double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) -integer, intent(in) :: lmax,kmax,n_kl(kmax,0:lmax) -integer, intent(in) :: n_a(3),n_b(3) -double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) - -! -! | _ _ _. | -! |_ (_) (_ (_| | -! - -double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm -double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big -double precision :: areal,freal,breal,t1,t2,int_prod_bessel -double precision :: arg - -integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot -integer :: l,k, nkl_max - -! _ -! |_) o _ _. ._ ._ _. -! |_) | (_| (_| | | (_| \/ -! _| / - -double precision, allocatable :: array_coefs_A(:,:) -double precision, allocatable :: array_coefs_B(:,:) - -double precision, allocatable :: array_R(:,:,:,:,:) -double precision, allocatable :: array_I_A(:,:,:,:,:) -double precision, allocatable :: array_I_B(:,:,:,:,:) - -double precision :: f1, f2, f3 - -if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then - Vpseudo=0.d0 - return -end if - -fourpi=4.d0*dacos(-1.d0) -ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) -bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) -arg= g_a*ac*ac + g_b*bc*bc - -if(arg.gt.-dlog(1.d-20))then - Vpseudo=0.d0 - return -endif - -freal=dexp(-arg) - -areal=2.d0*g_a*ac -breal=2.d0*g_b*bc -ntotA=n_a(1)+n_a(2)+n_a(3) -ntotB=n_b(1)+n_b(2)+n_b(3) -ntot=ntotA+ntotB - -nkl_max=4 - -allocate (array_coefs_A(0:ntot,3)) -allocate (array_coefs_B(0:ntot,3)) - -allocate (array_R(kmax,0:ntot+nkl_max,0:lmax,0:lmax+ntot,0:lmax+ntot)) - -allocate (array_I_A(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot)) - -allocate (array_I_B(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot)) - -if(ac.eq.0.d0.and.bc.eq.0.d0)then - - - accu=0.d0 - - do k=1,kmax - do l=0,lmax - ktot=ntot+n_kl(k,l) - if (v_kl(k,l) == 0.d0) cycle - do m=-l,l - prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) - if (prod == 0.d0) cycle - prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - if (prodp == 0.d0) cycle - accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) - enddo - enddo - enddo - - Vpseudo=accu*fourpi - -else if(ac.ne.0.d0.and.bc.ne.0.d0)then - - f=fourpi*fourpi - - theta_AC0=dacos( (a(3)-c(3))/ac ) - phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) - theta_BC0=dacos( (b(3)-c(3))/bc ) - phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) - - - do lambdap=0,lmax+ntotB - do lambda=0,lmax+ntotA - do l=0,lmax - do ktot=0,ntotA+ntotB+nkl_max - do k=1,kmax - array_R(k,ktot,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg) - enddo - enddo - enddo - enddo - enddo - - do k1=0,n_a(1) - array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1) - enddo - do k2=0,n_a(2) - array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2) - enddo - do k3=0,n_a(3) - array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3) - enddo - - do k1p=0,n_b(1) - array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p) - enddo - do k2p=0,n_b(2) - array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p) - enddo - do k3p=0,n_b(3) - array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p) - enddo - - accu=0.d0 - do l=0,lmax - do m=-l,l - - do k3=0,n_a(3) - do k2=0,n_a(2) - do k1=0,n_a(1) - do lambda=0,l+ntotA - do mu=-lambda,lambda - array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) - enddo - enddo - enddo - enddo - enddo - - do k3p=0,n_b(3) - do k2p=0,n_b(2) - do k1p=0,n_b(1) - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) - enddo - enddo - enddo - enddo - enddo - - do k3=0,n_a(3) - if (array_coefs_A(k3,3) == 0.d0) cycle - do k2=0,n_a(2) - if (array_coefs_A(k2,2) == 0.d0) cycle - do k1=0,n_a(1) - if (array_coefs_A(k1,1) == 0.d0) cycle - - do lambda=0,l+ntotA - do mu=-lambda,lambda - - prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*array_I_A(mu,lambda,k1,k2,k3) - if (prod == 0.d0) cycle - - do k3p=0,n_b(3) - do k2p=0,n_b(2) - do k1p=0,n_b(1) - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - - prodp=prod*ylm(lambdap,mup,theta_BC0,phi_BC0)* & - array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* & - array_I_B(mup,lambdap,k1p,k2p,k3p) - - if (prodp == 0.d0) cycle - do k=1,kmax - ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l) - accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap) - enddo - - enddo - enddo - enddo - - enddo - enddo - - enddo - enddo - enddo - - enddo - enddo - - enddo - enddo - - Vpseudo=f*accu - -else if(ac.eq.0.d0.and.bc.ne.0.d0)then - - f=fourpi**1.5d0 - theta_BC0=dacos( (b(3)-c(3))/bc ) - phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) - - areal=2.d0*g_a*ac - breal=2.d0*g_b*bc - freal=dexp(-g_a*ac**2-g_b*bc**2) - - do lambdap=0,lmax+ntotB - do l=0,lmax - do ktot=0,ntotA+ntotB+nkl_max - do k=1,kmax - array_R(k,ktot,l,0,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg) - enddo - enddo - enddo - enddo - - do k1p=0,n_b(1) - array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p) - enddo - do k2p=0,n_b(2) - array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p) - enddo - do k3p=0,n_b(3) - array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p) - enddo - - accu=0.d0 - do l=0,lmax - do m=-l,l - - do k3p=0,n_b(3) - do k2p=0,n_b(2) - do k1p=0,n_b(1) - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) - enddo - enddo - enddo - enddo - enddo - - prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) - - do k3p=0,n_b(3) - if (array_coefs_B(k3p,3) == 0.d0) cycle - do k2p=0,n_b(2) - if (array_coefs_B(k2p,2) == 0.d0) cycle - do k1p=0,n_b(1) - if (array_coefs_B(k1p,1) == 0.d0) cycle - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - - prodp=prod*array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(mup,lambdap,k1p,k2p,k3p) - - if (prodp == 0.d0) cycle - - do k=1,kmax - - ktot=ntotA+k1p+k2p+k3p+n_kl(k,l) - accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,0,lambdap) - - enddo - - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - Vpseudo=f*accu - -else if(ac.ne.0.d0.and.bc.eq.0.d0)then - - f=fourpi**1.5d0 - theta_AC0=dacos( (a(3)-c(3))/ac ) - phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) - - areal=2.d0*g_a*ac - breal=2.d0*g_b*bc - freal=dexp(-g_a*ac**2-g_b*bc**2) - - do lambda=0,lmax+ntotA - do l=0,lmax - do ktot=0,ntotA+ntotB+nkl_max - do k=1,kmax - array_R(k,ktot,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg) - enddo - enddo - enddo - enddo - - do k1=0,n_a(1) - array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1) - enddo - do k2=0,n_a(2) - array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2) - enddo - do k3=0,n_a(3) - array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3) - enddo - - accu=0.d0 - do l=0,lmax - do m=-l,l - - do k3=0,n_a(3) - do k2=0,n_a(2) - do k1=0,n_a(1) - do lambda=0,l+ntotA - do mu=-lambda,lambda - array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) - enddo - enddo - enddo - enddo - enddo - - do k3=0,n_a(3) - if (array_coefs_A(k3,3) == 0.d0) cycle - do k2=0,n_a(2) - if (array_coefs_A(k2,2) == 0.d0) cycle - do k1=0,n_a(1) - if (array_coefs_A(k1,1) == 0.d0) cycle - do lambda=0,l+ntotA - do mu=-lambda,lambda - - prod=array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(mu,lambda,k1,k2,k3) - if (prod == 0.d0) cycle - prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - - if (prodp == 0.d0) cycle - - do k=1,kmax - ktot=k1+k2+k3+ntotB+n_kl(k,l) - accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,0) - enddo - - enddo - enddo - enddo - enddo - enddo - - enddo - enddo - - Vpseudo=f*accu -endif - -! _ -! |_ o ._ _. | o _ _ -! | | | | (_| | | _> (/_ -! - deallocate (array_R, array_I_A, array_I_B) - deallocate (array_coefs_A, array_coefs_B) - return -end - -! _ -! | | -!__ __ _ __ ___ ___ _ _ __| | ___ _ __ _ _ _ __ ___ -!\ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ | '_ \| | | | '_ ` _ \ -! \ V / | |_) \__ \ __/ |_| | (_| | (_) | | | | | |_| | | | | | | -! \_/ | .__/|___/\___|\__,_|\__,_|\___/ |_| |_|\__,_|_| |_| |_| -! | | -! |_| - -double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) -implicit none - - -! ___ -! | ._ ._ _|_ -! _|_ | | |_) |_| |_ -! | -double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) -integer, intent(in) :: lmax,kmax,npts -integer, intent(in) :: n_a(3),n_b(3), n_kl(kmax,0:lmax) -double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) -double precision, intent(in) :: rmax - -! -! | _ _ _. | -! |_ (_) (_ (_| | -! - -integer :: l,m,k,kk -double precision ac(3),bc(3) -double precision dr,sum,rC -double precision overlap_orb_ylm_brute - -! _ -! / _. | _ | -! \_ (_| | (_ |_| | -! - -do l=1,3 - ac(l)=a(l)-c(l) - bc(l)=b(l)-c(l) -enddo - -dr=rmax/npts -sum=0.d0 -do l=0,lmax - do m=-l,l - do k=1,npts - rC=(k-1)*dr+dr/2.d0 - do kk=1,kmax - sum=sum+dr*v_kl(kk,l)*rC**(n_kl(kk,l)+2)*dexp(-dz_kl(kk,l)*rC**2) & - *overlap_orb_ylm_brute(npts,rC,n_a,ac,g_a,l,m) & - *overlap_orb_ylm_brute(npts,rC,n_b,bc,g_b,l,m) - enddo - enddo - enddo -enddo -Vpseudo_num=sum -return -end -!! Routine Vloc is a variation of formumla (66) -!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976) -!! without the projection operator -!! -!! Vloc= (4pi)**3/2* \sum_{k=1}^{klocmax} \sum_{l=0}^lmax \sum_{m=-l}^{l} -!!\sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} -!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3) -!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* -!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} -!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p) -!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* -!!\sum_{l=0}^lmax \sum_{m=-l}^{l} - -!! bigI(0,0,l,m,k1+k1p,k2+k2p,k3+k3p)*Y_{l m}(D_unit) -!! *v_k(k)* bigR(lambda,k1+k2+k3+k1p+k2p+k3p+n_k(k),g_a,g_b,AC,BC,dz_k(k)) -!! -!! nA=nAx+nAy+nAz -!! nB=nBx+nBy+nBz -!! D=(g_a AC+g_b BC) -!! D_unit= vect(D)/D -!! AC_x= A_x-C_x, etc. -!! BC=|B-C| -!! AC_unit= vect(AC)/AC -!! BC_unit= vect(BC)/BCA -!! -!! bigR(lambda,g_a,g_b,g_k,AC,BC) -!! = exp(-g_a* AC**2 -g_b* BC**2)* -!! I_loc= \int dx x**l *exp(-gam*x**2) M_n(ax) l=ktot+2 gam=g_a+g_b+dz_k(k) a=dreal n=l -!! M_n(x) modified spherical bessel function - - -double precision function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) -implicit none -integer klocmax -double precision v_k(klocmax),dz_k(klocmax),crochet,bigA -integer n_k(klocmax) -double precision a(3),g_a,b(3),g_b,c(3),d(3) -integer n_a(3),n_b(3),ntotA,ntotB,ntot,m -integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p -double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0,coef -double precision,allocatable :: array_R_loc(:,:,:) -double precision,allocatable :: array_coefs(:,:,:,:,:,:) -double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg - - - fourpi=4.d0*dacos(-1.d0) - f=fourpi**1.5d0 - ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) - bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) - arg=g_a*ac**2+g_b*bc**2 - if(arg.gt.-dlog(1.d-20))then - Vloc=0.d0 - return - endif - - ntotA=n_a(1)+n_a(2)+n_a(3) - ntotB=n_b(1)+n_b(2)+n_b(3) - ntot=ntotA+ntotB - - if(ac.eq.0.d0.and.bc.eq.0.d0)then - accu=0.d0 - - do k=1,klocmax - accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) - enddo - Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) - !bigI frequently is null - return - endif - - freal=dexp(-g_a*ac**2-g_b*bc**2) - - d2 = 0.d0 - do i=1,3 - d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) - d2=d2+d(i)*d(i) - enddo - d2=dsqrt(d2) - dreal=2.d0*d2 - - - allocate (array_R_loc(-2:ntot+klocmax,klocmax,0:ntot)) - allocate (array_coefs(0:ntot,0:ntot,0:ntot,0:ntot,0:ntot,0:ntot)) - - do ktot=-2,ntotA+ntotB+klocmax - do l=0,ntot - do k=1,klocmax - array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) - enddo - enddo - enddo - - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)& - *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)& - *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)& - *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) - enddo - enddo - enddo - enddo - enddo - enddo - - - accu=0.d0 - if(d2 == 0.d0)then - l=0 - m=0 - coef=1.d0/dsqrt(4.d0*dacos(-1.d0)) - do k=1,klocmax - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & - *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) - ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) - accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - else - theta_DC0=dacos(d(3)/d2) - phi_DC0=datan2(d(2)/d2,d(1)/d2) - - do k=1,klocmax - if (v_k(k) == 0.d0) cycle - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - if (array_coefs(k1,k2,k3,k1p,k2p,k3p) == 0.d0) cycle - do l=0,ntot - do m=-l,l - coef=ylm(l,m,theta_DC0,phi_DC0) - if (coef == 0.d0) cycle - ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) - if (array_R_loc(ktot,k,l) == 0.d0) cycle - prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & - *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) - accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - enddo - enddo - endif - Vloc=f*accu - - deallocate (array_R_loc) - deallocate (array_coefs) -end - -double precision function bigA(i,j,k) -implicit none -integer i,j,k -double precision fourpi,dble_fact -fourpi=4.d0*dacos(-1.d0) -bigA=0.d0 -if(mod(i,2).eq.1)return -if(mod(j,2).eq.1)return -if(mod(k,2).eq.1)return -bigA=fourpi*dble_fact(i-1)*dble_fact(j-1)*dble_fact(k-1)/dble_fact(i+j+k+1) -end -!! -!! I_{lambda,mu,l,m}^{k1,k2,k3} = /int dOmega Y_{lambda mu} xchap^k1 ychap^k2 zchap^k3 Y_{lm} -!! - -double precision function bigI(lambda,mu,l,m,k1,k2,k3) -implicit none -integer lambda,mu,l,m,k1,k2,k3 -integer k,i,kp,ip -double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm -double precision sgn, sgnp -pi=dacos(-1.d0) - -bigI=0.d0 -if(mu.gt.0.and.m.gt.0)then -sum=0.d0 -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2== 0.d0) return -sgn = 1.d0 -do k=0,mu/2 - do i=0,lambda-mu - if (coef_pm(lambda,i+mu) == 0.d0) cycle - sgnp = 1.d0 - do kp=0,m/2 - do ip=0,l-m - cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo - enddo - sgn = -sgn -enddo -bigI=sum -return -endif - -if(mu.eq.0.and.m.eq.0)then -factor1=dsqrt((2*lambda+1)/(4.d0*pi)) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)/(4.d0*pi)) -if (factor2== 0.d0) return -sum=0.d0 -do i=0,lambda - do ip=0,l - cylm=factor1*coef_pm(lambda,i) - if (cylm == 0.d0) cycle - cylmp=factor2*coef_pm(l,ip) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(k1,k2,i+ip+k3) - enddo -enddo -bigI=sum -return -endif - -if(mu.eq.0.and.m.gt.0)then -factor1=dsqrt((2*lambda+1)/(4.d0*pi)) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2== 0.d0) return -sum=0.d0 -do i=0,lambda - sgnp = 1.d0 - do kp=0,m/2 - do ip=0,l-m - cylm=factor1*coef_pm(lambda,i) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo -enddo -bigI=sum -return -endif - -if(mu.gt.0.and.m.eq.0)then -sum=0.d0 -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)/(4.d0*pi)) -if (factor2== 0.d0) return -sgn = 1.d0 -do k=0,mu/2 - do i=0,lambda-mu - if (coef_pm(lambda,i+mu) == 0.d0) cycle - do ip=0,l - cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=factor2*coef_pm(l,ip) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3) - enddo - enddo - sgn = -sgn -enddo -bigI=sum -return -endif - -if(mu.lt.0.and.m.lt.0)then -mu=-mu -m=-m -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2== 0.d0) return -sum=0.d0 -sgn = 1.d0 -do k=0,(mu-1)/2 - do i=0,lambda-mu - if (coef_pm(lambda,i+mu) == 0.d0) cycle - sgnp = 1.d0 - do kp=0,(m-1)/2 - do ip=0,l-m - if (coef_pm(l,ip+m) == 0.d0) cycle - cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo - enddo - sgn = -sgn -enddo -mu=-mu -m=-m -bigI=sum -return -endif - -if(mu.eq.0.and.m.lt.0)then -m=-m -factor1=dsqrt((2*lambda+1)/(4.d0*pi)) -if (factor1 == 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2 == 0.d0) return -sum=0.d0 -do i=0,lambda - sgnp = 1.d0 - do kp=0,(m-1)/2 - do ip=0,l-m - cylm=factor1*coef_pm(lambda,i) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo -enddo -m=-m -bigI=sum -return -endif - -if(mu.lt.0.and.m.eq.0)then -sum=0.d0 -mu=-mu -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)/(4.d0*pi)) -if (factor2== 0.d0) return -sgn = 1.d0 -do k=0,(mu-1)/2 - do i=0,lambda-mu - do ip=0,l - cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=factor2*coef_pm(l,ip) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3) - enddo - enddo - sgn = -sgn -enddo -mu=-mu -bigI=sum -return -endif - -if(mu.gt.0.and.m.lt.0)then -sum=0.d0 -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2== 0.d0) return -m=-m -sgn=1.d0 -do k=0,mu/2 - do i=0,lambda-mu - if (coef_pm(lambda,i+mu) == 0.d0) cycle - sgnp=1.d0 - do kp=0,(m-1)/2 - do ip=0,l-m - if (coef_pm(l,ip+m) == 0.d0) cycle - cylm =sgn *factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo - enddo - sgn = -sgn -enddo -m=-m -bigI=sum -return -endif - -if(mu.lt.0.and.m.gt.0)then -mu=-mu -factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) -if (factor1== 0.d0) return -factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) -if (factor2== 0.d0) return -sum=0.d0 -sgn = 1.d0 -do k=0,(mu-1)/2 - do i=0,lambda-mu - if (coef_pm(lambda,i+mu) == 0.d0) cycle - sgnp = 1.d0 - do kp=0,m/2 - do ip=0,l-m - if (coef_pm(l,ip+m) == 0.d0) cycle - cylm=sgn*factor1 *binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) - if (cylm == 0.d0) cycle - cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) - if (cylmp == 0.d0) cycle - sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3) - enddo - sgnp = -sgnp - enddo - enddo - sgn = -sgn -enddo -bigI=sum -mu=-mu -return -endif - -stop 'pb in bigI!' -end - -double precision function crochet(n,g) -implicit none -integer n -double precision g,dble_fact,expo -double precision, parameter :: sq_pi_ov_2=dsqrt(dacos(-1.d0)*0.5d0) -expo=0.5d0*dfloat(n+1) -crochet=dble_fact(n-1)/(g+g)**expo -if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2 -end - -!! -!! overlap= = /int dOmega Ylm (x-center_x)**nx*(y-center_y)**nx*(z-center)**nx -!! *exp(-g*(r-center)**2) -!! -double precision function overlap_orb_ylm_brute(npts,r,npower_orb,center_orb,g_orb,l,m) -implicit none -integer npower_orb(3),l,m,i,j,npts -double precision u,g_orb,du,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) -double precision x_orb,y_orb,z_orb,twopi,r -twopi=2.d0*dacos(-1.d0) -du=2.d0/npts -dphi=twopi/npts -overlap_orb_ylm_brute=0.d0 -do i=1,npts - u=-1.d0+du*(i-1)+du/2.d0 - sintheta=dsqrt(1.d0-u**2) - do j=1,npts - phi=dphi*(j-1)+dphi/2.d0 - x_orb=r*dcos(phi)*sintheta - y_orb=r*dsin(phi)*sintheta - z_orb=r*u - term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,u,phi) - overlap_orb_ylm_brute= overlap_orb_ylm_brute+term*du*dphi - enddo -enddo -end - -double precision function overlap_orb_ylm_grid(nptsgrid,r_orb,npower_orb,center_orb,g_orb,l,m) -implicit none -!! PSEUDOS -integer nptsgridmax,nptsgrid -double precision coefs_pseudo,ptsgrid -parameter(nptsgridmax=50) -common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) -!!!!! -integer npower_orb(3),l,m,i -double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) -double precision x_orb,y_orb,z_orb,twopi,pi,cosphi,sinphi,xbid -pi=dacos(-1.d0) -twopi=2.d0*pi -overlap_orb_ylm_grid=0.d0 -do i=1,nptsgrid - x_orb=r_orb*ptsgrid(i,1) - y_orb=r_orb*ptsgrid(i,2) - z_orb=r_orb*ptsgrid(i,3) - x=ptsgrid(i,3) - phi=datan2(ptsgrid(i,2),ptsgrid(i,1)) - term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,x,phi) - overlap_orb_ylm_grid= overlap_orb_ylm_grid+coefs_pseudo(i)*term -enddo -overlap_orb_ylm_grid=2.d0*twopi*overlap_orb_ylm_grid -end - -! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) exp(i m phi) -! l=0,1,2,.... -! m=0,1,...,l -! Here: -! n=l (n=0,1,...) -! m=0,1,...,n -! x=cos(theta) 0 < x < 1 -! -! -! This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n - -! Exemples (see 'Associated Legendre Polynomilas wikipedia') -! P_{0}^{0}(x)=1 -! P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) -! P_{1}^{0}(x)=x -! P_{1}^{1}(x)=-(1-x^2)^{1/2} -! P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) -! P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) -! P_{2}^{0}(x)=1/2 (3x^{2}-1) -! P_{2}^{1}(x)=-3x(1-x^2)^{1/2} -! P_{2}^{2}(x)=3(1-x^2) - - - SUBROUTINE LPMN(MM,M,N,X,PM) -! -! Here N = LMAX -! Here M= MMAX (we take M=LMAX in input) -! -! ===================================================== -! Purpose: Compute the associated Legendre functions Pmn(x) -! Input : x --- Argument of Pmn(x) -! m --- Order of Pmn(x), m = 0,1,2,...,n -! n --- Degree of Pmn(x), n = 0,1,2,...,N -! mm --- Physical dimension of PM -! Output: PM(m,n) --- Pmn(x) -! ===================================================== -! - implicit none -! IMPLICIT DOUBLE PRECISION (P,X) - integer :: MM, N, I, J, M - double precision :: PM(0:MM,0:(N+1)), X, XQ, XS - DOUBLE PRECISION, SAVE :: INVERSE(100) = 0.D0 - DOUBLE PRECISION :: LS, II, JJ - IF (INVERSE(1) == 0.d0) THEN - DO I=1,100 - INVERSE(I) = 1.D0/DBLE(I) - ENDDO - ENDIF - DO I=0,N - DO J=0,M - PM(J,I)=0.0D0 - ENDDO - ENDDO - PM(0,0)=1.0D0 - IF (DABS(X).EQ.1.0D0) THEN - DO I=1,N - PM(0,I)=X**I - ENDDO - RETURN - ENDIF - LS=1.D0 - IF (DABS(X).GT.1.0D0) LS=-1.D0 - XQ=DSQRT(LS*(1.0D0-X*X)) - XS=LS*(1.0D0-X*X) - II = 1.D0 - DO I=1,M - PM(I,I)=-LS*II*XQ*PM(I-1,I-1) - II = II+2.D0 - ENDDO - II = 1.D0 - DO I=0,M - PM(I,I+1)=II*X*PM(I,I) - II = II+2.D0 - ENDDO - - II = 0.D0 - DO I=0,M - JJ = II+2.D0 - DO J=I+2,N - PM(I,J)=((2.0D0*JJ-1.0D0)*X*PM(I,J-1)- (II+JJ-1.0D0)*PM(I,J-2))*INVERSE(J-I) - JJ = JJ+1.D0 - ENDDO - II = II+1.D0 - ENDDO - END - - -! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 -! P_l^|m|(cos(theta)) exp(i m phi) - - subroutine erreur(x,n,rmoy,error) - implicit none - integer :: i, n - double precision :: x(n), rn, rn1, error, rmoy -! calcul de la moyenne - rmoy=0.d0 - do i=1,n - rmoy=rmoy+x(i) - enddo - rmoy=rmoy/dfloat(n) -! calcul de l'erreur - error=0.d0 - do i=1,n - error=error+(x(i)-rmoy)**2 - enddo - if(n.gt.1)then - rn=dfloat(n) - rn1=dfloat(n-1) - error=dsqrt(error)/dsqrt(rn*rn1) - else - write(2,*)'Seulement un block Erreur nondefinie' - error=0.d0 - endif - end - - subroutine initpseudos(nptsgrid) - implicit none - integer nptsgridmax,nptsgrid,ik - double precision coefs_pseudo,ptsgrid - double precision p,q,r,s - parameter(nptsgridmax=50) - common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) - - p=1.d0/dsqrt(2.d0) - q=1.d0/dsqrt(3.d0) - r=1.d0/dsqrt(11.d0) - s=3.d0/dsqrt(11.d0) - - if(nptsgrid.eq.4)then - - ptsgrid(1,1)=q - ptsgrid(1,2)=q - ptsgrid(1,3)=q - - ptsgrid(2,1)=q - ptsgrid(2,2)=-q - ptsgrid(2,3)=-q - - ptsgrid(3,1)=-q - ptsgrid(3,2)=q - ptsgrid(3,3)=-q - - ptsgrid(4,1)=-q - ptsgrid(4,2)=-q - ptsgrid(4,3)=q - - do ik=1,4 - coefs_pseudo(ik)=1.d0/4.d0 - enddo - return - endif - - ptsgrid(1,1)=1.d0 - ptsgrid(1,2)=0.d0 - ptsgrid(1,3)=0.d0 - - ptsgrid(2,1)=-1.d0 - ptsgrid(2,2)=0.d0 - ptsgrid(2,3)=0.d0 - - ptsgrid(3,1)=0.d0 - ptsgrid(3,2)=1.d0 - ptsgrid(3,3)=0.d0 - - ptsgrid(4,1)=0.d0 - ptsgrid(4,2)=-1.d0 - ptsgrid(4,3)=0.d0 - - ptsgrid(5,1)=0.d0 - ptsgrid(5,2)=0.d0 - ptsgrid(5,3)=1.d0 - - ptsgrid(6,1)=0.d0 - ptsgrid(6,2)=0.d0 - ptsgrid(6,3)=-1.d0 - - do ik=1,6 - coefs_pseudo(ik)=1.d0/6.d0 - enddo - - if(nptsgrid.eq.6)return - - ptsgrid(7,1)=p - ptsgrid(7,2)=p - ptsgrid(7,3)=0.d0 - - ptsgrid(8,1)=p - ptsgrid(8,2)=-p - ptsgrid(8,3)=0.d0 - - ptsgrid(9,1)=-p - ptsgrid(9,2)=p - ptsgrid(9,3)=0.d0 - - ptsgrid(10,1)=-p - ptsgrid(10,2)=-p - ptsgrid(10,3)=0.d0 - - ptsgrid(11,1)=p - ptsgrid(11,2)=0.d0 - ptsgrid(11,3)=p - - ptsgrid(12,1)=p - ptsgrid(12,2)=0.d0 - ptsgrid(12,3)=-p - - ptsgrid(13,1)=-p - ptsgrid(13,2)=0.d0 - ptsgrid(13,3)=p - - ptsgrid(14,1)=-p - ptsgrid(14,2)=0.d0 - ptsgrid(14,3)=-p - - ptsgrid(15,1)=0.d0 - ptsgrid(15,2)=p - ptsgrid(15,3)=p - - ptsgrid(16,1)=0.d0 - ptsgrid(16,2)=p - ptsgrid(16,3)=-p - - ptsgrid(17,1)=0.d0 - ptsgrid(17,2)=-p - ptsgrid(17,3)=p - - ptsgrid(18,1)=0.d0 - ptsgrid(18,2)=-p - ptsgrid(18,3)=-p - - do ik=1,6 - coefs_pseudo(ik)=1.d0/30.d0 - enddo - do ik=7,18 - coefs_pseudo(ik)=1.d0/15.d0 - enddo - - if(nptsgrid.eq.18)return - - ptsgrid(19,1)=q - ptsgrid(19,2)=q - ptsgrid(19,3)=q - - ptsgrid(20,1)=-q - ptsgrid(20,2)=q - ptsgrid(20,3)=q - - ptsgrid(21,1)=q - ptsgrid(21,2)=-q - ptsgrid(21,3)=q - - ptsgrid(22,1)=q - ptsgrid(22,2)=q - ptsgrid(22,3)=-q - - ptsgrid(23,1)=-q - ptsgrid(23,2)=-q - ptsgrid(23,3)=q - - ptsgrid(24,1)=-q - ptsgrid(24,2)=q - ptsgrid(24,3)=-q - - ptsgrid(25,1)=q - ptsgrid(25,2)=-q - ptsgrid(25,3)=-q - - ptsgrid(26,1)=-q - ptsgrid(26,2)=-q - ptsgrid(26,3)=-q - - do ik=1,6 - coefs_pseudo(ik)=1.d0/21.d0 - enddo - do ik=7,18 - coefs_pseudo(ik)=4.d0/105.d0 - enddo - do ik=19,26 - coefs_pseudo(ik)=27.d0/840.d0 - enddo - - if(nptsgrid.eq.26)return - - ptsgrid(27,1)=r - ptsgrid(27,2)=r - ptsgrid(27,3)=s - - ptsgrid(28,1)=r - ptsgrid(28,2)=-r - ptsgrid(28,3)=s - - ptsgrid(29,1)=-r - ptsgrid(29,2)=r - ptsgrid(29,3)=s - - ptsgrid(30,1)=-r - ptsgrid(30,2)=-r - ptsgrid(30,3)=s - - ptsgrid(31,1)=r - ptsgrid(31,2)=r - ptsgrid(31,3)=-s - - ptsgrid(32,1)=r - ptsgrid(32,2)=-r - ptsgrid(32,3)=-s - - ptsgrid(33,1)=-r - ptsgrid(33,2)=r - ptsgrid(33,3)=-s - - ptsgrid(34,1)=-r - ptsgrid(34,2)=-r - ptsgrid(34,3)=-s - - ptsgrid(35,1)=r - ptsgrid(35,2)=s - ptsgrid(35,3)=r - - ptsgrid(36,1)=-r - ptsgrid(36,2)=s - ptsgrid(36,3)=r - - ptsgrid(37,1)=r - ptsgrid(37,2)=s - ptsgrid(37,3)=-r - - ptsgrid(38,1)=-r - ptsgrid(38,2)=s - ptsgrid(38,3)=-r - - ptsgrid(39,1)=r - ptsgrid(39,2)=-s - ptsgrid(39,3)=r - - ptsgrid(40,1)=r - ptsgrid(40,2)=-s - ptsgrid(40,3)=-r - - ptsgrid(41,1)=-r - ptsgrid(41,2)=-s - ptsgrid(41,3)=r - - ptsgrid(42,1)=-r - ptsgrid(42,2)=-s - ptsgrid(42,3)=-r - - ptsgrid(43,1)=s - ptsgrid(43,2)=r - ptsgrid(43,3)=r - - ptsgrid(44,1)=s - ptsgrid(44,2)=r - ptsgrid(44,3)=-r - - ptsgrid(45,1)=s - ptsgrid(45,2)=-r - ptsgrid(45,3)=r - - ptsgrid(46,1)=s - ptsgrid(46,2)=-r - ptsgrid(46,3)=-r - - ptsgrid(47,1)=-s - ptsgrid(47,2)=r - ptsgrid(47,3)=r - - ptsgrid(48,1)=-s - ptsgrid(48,2)=r - ptsgrid(48,3)=-r - - ptsgrid(49,1)=-s - ptsgrid(49,2)=-r - ptsgrid(49,3)=r - - ptsgrid(50,1)=-s - ptsgrid(50,2)=-r - ptsgrid(50,3)=-r - - do ik=1,6 - coefs_pseudo(ik)=4.d0/315.d0 - enddo - do ik=7,18 - coefs_pseudo(ik)=64.d0/2835.d0 - enddo - do ik=19,26 - coefs_pseudo(ik)=27.d0/1280.d0 - enddo - do ik=27,50 - coefs_pseudo(ik)=14641.d0/725760.d0 - enddo - - if(nptsgrid.eq.50)return - - write(*,*)'Grid for pseudos not available!' - write(*,*)'N=4-6-18-26-50 only!' - stop - end - -!! -!! R_{lambda,lamba',N}= exp(-ga_a AC**2 -g_b BC**2) /int_{0}{+infty} r**(2+n) exp(-(g_a+g_b+g_k)r**2) -!! * M_{lambda}( 2g_a ac r) M_{lambda'}(2g_b bc r) -!! - double precision function bigR(lambda,lambdap,n,g_a,g_b,ac,bc,g_k) - implicit none - integer lambda,lambdap,n,npts,i - double precision g_a,g_b,ac,bc,g_k,arg,factor,delta1,delta2,cc,rmax,dr,sum,x1,x2,r - double precision bessel_mod - arg=g_a*ac**2+g_b*bc**2 - factor=dexp(-arg) - delta1=2.d0*g_a*ac - delta2=2.d0*g_b*bc - cc=g_a+g_b+g_k - if(cc.eq.0.d0)stop 'pb. in bigR' - rmax=dsqrt(-dlog(10.d-20)/cc) - npts=500 - dr=rmax/npts - sum=0.d0 - do i=1,npts - r=(i-1)*dr - x1=delta1*r - x2=delta2*r - sum=sum+dr*r**(n+2)*dexp(-cc*r*r)*bessel_mod(x1,lambda)*bessel_mod(x2,lambdap) - enddo - bigR=sum*factor - end - - double precision function bessel_mod(x,n) - implicit none - integer n - double precision x,bessel_mod_exp,bessel_mod_recur - if(x.le.0.8d0)then - bessel_mod=bessel_mod_exp(n,x) - else - bessel_mod=bessel_mod_recur(n,x) - endif - end - - recursive function bessel_mod_recur(n,x) result(a) - implicit none - integer n - double precision x,a,bessel_mod_exp - if(x.le.0.8d0)then - a=bessel_mod_exp(n,x) - return - endif - if(n.eq.0)a=dsinh(x)/x - if(n.eq.1)a=(x*dcosh(x)-dsinh(x))/(x*x) - if(n.ge.2)a=bessel_mod_recur(n-2,x)-(n+n-1)/x*bessel_mod_recur(n-1,x) - end - - double precision function bessel_mod_exp(n,x) - implicit none - integer n,k - double precision x,coef,accu,fact,dble_fact - accu=0.d0 - do k=0,10 - coef=1.d0/(fact(k)*dble_fact(2*(n+k)+1)) - accu=accu+(0.5d0*x*x)**k*coef - enddo - bessel_mod_exp=x**n*accu - end - - - - -!c Computation of real spherical harmonics Ylm(theta,phi) -!c -!c l=0,1,.... -!c m=-l,l -!c -!c m>0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) cos(m phi) -!c m=0: Y_l0 = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^0 (cos(theta)) -!c m<0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) sin(|m|phi) - -!Examples(wikipedia http://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics) - -! l = 0 - -! Y_00 = \sqrt{1/(4pi)} - -! l = 1 - -! Y_1-1= \sqrt{3/(4pi)} y/r -! Y_10 = \sqrt{3/(4pi)} z/r -! Y_11 = \sqrt{3/(4pi)} x/r -! -! l = 2 -! -! Y_2,-2= 1/2 \sqrt{15/pi} xy/r^2 -! Y_2,-1= 1/2 \sqrt{15/pi} yz/r^2 -! Y_20 = 1/4 \sqrt{15/pi} (-x^2-y^2 +2z^2)/r^2 -! Y_21 = 1/2 \sqrt{15/pi} zx/r^2 -! Y_22 = 1/4 \sqrt{15/pi} (x^2-y^2)/r^2 -! -!c -double precision function ylm(l,m,theta,phi) -implicit none -integer l,m,i -double precision theta,phi,pm,factor,twopi,x,fact,sign -DIMENSION PM(0:100,0:100) -twopi=2.d0*dacos(-1.d0) -x=dcos(theta) -if (iand(m,1) == 1) then - sign=-1.d0 -else - sign=1.d0 -endif -CALL LPMN(100,l,l,X,PM) -if (m > 0) then - factor=dsqrt((l+l+1)*fact(l-m) /(twopi*fact(l+m)) ) -! factor = dble(l+m) -! do i=-m,m-1 -! factor = factor * (factor - 1.d0) -! enddo -! factor=dsqrt(dble(l+l+1)/(twopi*factor) ) - ylm=sign*factor*pm(m,l)*dcos(dfloat(m)*phi) -else if (m == 0) then - factor=dsqrt( 0.5d0*(l+l+1) /twopi ) - ylm=factor*pm(m,l) -else if (m < 0) then - factor=dsqrt( (l+l+1)*fact(l+m) /(twopi*fact(l-m)) ) -! factor = dble(l-m) -! do i=m,-m-1 -! factor = factor * (factor - 1.d0) -! enddo -! factor=dsqrt(dble(l+l+1)/(twopi*factor) ) - ylm=sign*factor*pm(-m,l)*dsin(dfloat(-m)*phi) -endif -end - -!c Explicit representation of Legendre polynomials P_n(x) -!! -!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k -!! -!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) -!! coef_pm(n,k) is the k_th coefficient of P_n(x) -double precision function coef_pm(n,k) -implicit none -integer n,k -double precision arg,binom_func,binom_gen -if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined' -if(n.eq.0.and.k.eq.0)then -coef_pm=1.d0 -return -endif -arg=0.5d0*dfloat(n+k-1) -coef_pm=2.d0**n*binom_func(n,k)*binom_gen(arg,n) -end - -!! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k -!! xchap=x/r etc. -!c m>0: Y_lm = sqrt(2)*factor* P_l^|m|(cos(theta)) cos(m phi) -!c m=0: Y_l0 = factor* P_l^0 (cos(theta)) -!c m<0: Y_lm = sqrt(2) factor* P_l^|m|(cos(theta)) sin(|m|phi) -!c factor= ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 - -!! P_l^m (x) = (-1)**m (1-x**2)^m/2 d^m/dx^m P_l(x) m >0 or 0 -!! the series expansion of P_m (x) is known -!! -!! sin(theta)**m cos(mphi) = \sum_0^[m/2] binom(m,2k) x^(m-2k) y^2k (-1)**k (easy to proove with -!! Moivre formula) -!! (here x = xchap...) -!! -!! Ylm m> 0 = \sum_{k=0}^[m/2] \sum_{i=0}^(l-m) c_ki x^(m-2k) y^2k z^i -!! -!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k)*(m+i)!/i!*coef_pm(l,i+m) -!! -!! Ylm m< 0 = \sum_{k=0}^[(m-1)/2] \sum_{i=0}^(l-m) c_ki x^(m-(2k+1)) y^(2k+1) z^i -!! -!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k+1)*(m+i)!/i!*coef_pm(l,i+m) - - -double precision function ylm_bis(l,m,theta,phi) -implicit none -integer l,m,k,i -double precision x,y,z,theta,phi,sum,factor,pi,binom_func,fact,coef_pm,cylm -pi=dacos(-1.d0) -x=dsin(theta)*dcos(phi) -y=dsin(theta)*dsin(phi) -z=dcos(theta) -factor=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) -if(m.gt.0)then -sum=0.d0 -do k=0,m/2 - do i=0,l-m - cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m) - sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i - enddo -enddo -ylm_bis=sum -return -endif -if(m.eq.0)then -sum=0.d0 -do i=0,l - sum=sum+factor*coef_pm(l,i)*z**i -enddo -ylm_bis=sum -return -endif -if(m.lt.0)then -m=-m -sum=0.d0 -do k=0,(m-1)/2 - do i=0,l-m - cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m) - sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i - enddo -enddo -ylm_bis=sum -m=-m -return -endif -end - -!c -!c Computation of associated Legendre Polynomials PM(m,n) for n=0,...,N -!c Here: -!c n=l (n=0,1,...) -!c m=0,1,...,n -!c x=cos(theta) 0 < x < 1 -!c -!c This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n -!c Exemples (see 'Associated Legendre Polynomilas wikipedia') -!c P_{0}^{0}(x)=1 -!c P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) -!c P_{1}^{0}(x)=x -!c P_{1}^{1}(x)=-(1-x^2)^{1/2} -!c P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) -!c P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) -!c P_{2}^{0}(x)=1/2 (3x^{2}-1) -!c P_{2}^{1}(x)=-3x(1-x^2)^{1/2} -!c P_{2}^{2}(x)=3(1-x^2) -!c -!c Explicit representation: -!! -!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k -!! -!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) - -double precision function binom_gen(alpha,n) - implicit none - integer :: n,k - double precision :: fact,alpha,prod, factn_inv - prod=1.d0 - factn_inv = 1.d0/fact(n) - do k=0,n-1 - prod=prod*(alpha-k) - binom_gen = prod*factn_inv - enddo -end - - -double precision function coef_nk(n,k) - implicit none - integer n,k - - double precision gam,dble_fact,fact - - if (k<0) stop 'pseudopot.f90 : coef_nk' - if (k>63) stop 'pseudopot.f90 : coef_nk' - gam=dble_fact(n+n+k+k+1) -! coef_nk=1.d0/(2.d0**k*fact(k)*gam) - coef_nk=1.d0/(dble(ibset(0_8,k))*fact(k)*gam) - - return - -end - -!! Calculation of -!! -!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) M_m(bx) -!! -!! M_n(x) modified spherical bessel function -!! - -double precision function int_prod_bessel(l,gam,n,m,a,b,arg) - - implicit none - integer n,k,m,q,l,kcp - double precision gam,dble_fact,fact,pi,a,b - double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg - - integer :: n_1,m_1,nlm - double precision :: term_A, term_B, term_rap, expo - double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square - double precision :: int_prod_bessel_loc - double precision :: inverses(0:300) - double precision :: two_qkmp1, qk, mk, nk - - logical done - - u=(a+b)*0.5d0/dsqrt(gam) - freal=dexp(-arg) - - if(a.eq.0.d0.and.b.eq.0.d0)then - if(n.ne.0.or.m.ne.0)then - int_prod_bessel=0.d0 - return - endif - - int_prod_bessel=crochet(l,gam)*freal - return - endif - - if(u.gt.6.d0)then - int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg) - return - endif - - if(a.ne.0.d0.and.b.ne.0.d0)then - - q=0 - intold=-1.d0 - int=0.d0 - done=.false. - - n_1 = n+n+1 - m_1 = m+m+1 - nlm = n+m+l - pi=dacos(-1.d0) - a_over_b_square = (a*a)/(b*b) - - ! First term of the sequence - - term_a =dble_fact(nlm-1) / (dble_fact(n_1)*dble_fact(m_1)) - expo=0.5d0*dfloat(nlm+1) - term_rap = term_a / (2.d0*gam)**expo - - s_0_0=term_rap*a**(n)*b**(m) - if(mod(nlm,2).eq.0)s_0_0=s_0_0*dsqrt(pi*.5d0) - - ! Initialize the first recurrence term for the q loop - s_q_0 = s_0_0 - - - mk = dble(m) - ! Loop over q for the convergence of the sequence - do while (.not.done) - - ! Init - s_q_k=s_q_0 - sum=s_q_0 - - if (q>300) then - stop 'pseudopot.f90 : q > 300' - endif - - qk = dble(q) - two_qkmp1 = 2.d0*(qk+mk)+1.d0 - do k=0,q-1 - if (s_q_k < 1.d-32) then - s_q_k = 0.d0 - exit - endif - s_q_k = two_qkmp1*qk*inverses(k)*s_q_k - sum=sum+s_q_k - two_qkmp1 = two_qkmp1-2.d0 - qk = qk-1.d0 - enddo - inverses(q) = a_over_b_square/(dble(q+n+q+n+3) * dble(q+1)) - - int=int+sum - - if(dabs(int-intold).lt.1d-15)then - done=.true. - else - - !Compute the s_q+1_0 - s_q_0=s_q_0*(q+q+nlm+1)*b*b/(dble(8*(m+q)+12)*(q+1)*gam) - - if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi*.5d0) - ! Increment q - q=q+1 - intold=int - endif - - enddo - - int_prod_bessel=int*freal - - return - endif - - if(a.eq.0.d0.and.b.ne.0.d0)then - - int = int_prod_bessel_loc(l,gam,m,b) - int_prod_bessel=int*freal - return - endif - - if(a.ne.0.d0.and.b.eq.0.d0)then - - int = int_prod_bessel_loc(l,gam,n,a) - int_prod_bessel=int*freal - return - - endif - - stop 'pb in int_prod_bessel!!' -end - -double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg) - implicit none - integer n,m,i,npts,l - double precision gam,a,b - double precision sum,x,bessel_mod,u,factor,arg - double precision xq(100),wq(100) - - u=(a+b)/(2.d0*dsqrt(gam)) - factor=dexp(u*u-arg)/dsqrt(gam) - -xq(1)= 5.38748089001123 -xq(2)= 4.60368244955074 -xq(3)= 3.94476404011563 -xq(4)= 3.34785456738322 -xq(5)= 2.78880605842813 -xq(6)= 2.25497400208928 -xq(7)= 1.73853771211659 -xq(8)= 1.23407621539532 -xq(9)= 0.737473728545394 -xq(10)= 0.245340708300901 -xq(11)=-0.245340708300901 -xq(12)=-0.737473728545394 -xq(13)=-1.23407621539532 -xq(14)=-1.73853771211659 -xq(15)=-2.25497400208928 -xq(16)=-2.78880605842813 -xq(17)=-3.34785456738322 -xq(18)=-3.94476404011563 -xq(19)=-4.60368244955074 -xq(20)=-5.38748089001123 -wq(1)= 2.229393645534151E-013 -wq(2)= 4.399340992273176E-010 -wq(3)= 1.086069370769280E-007 -wq(4)= 7.802556478532063E-006 -wq(5)= 2.283386360163528E-004 -wq(6)= 3.243773342237853E-003 -wq(7)= 2.481052088746362E-002 -wq(8)= 0.109017206020022 -wq(9)= 0.286675505362834 -wq(10)= 0.462243669600610 -wq(11)= 0.462243669600610 -wq(12)= 0.286675505362834 -wq(13)= 0.109017206020022 -wq(14)= 2.481052088746362E-002 -wq(15)= 3.243773342237853E-003 -wq(16)= 2.283386360163528E-004 -wq(17)= 7.802556478532063E-006 -wq(18)= 1.086069370769280E-007 -wq(19)= 4.399340992273176E-010 -wq(20)= 2.229393645534151E-013 - - npts=20 -! call gauher(xq,wq,npts) - - sum=0.d0 - do i=1,npts - x=(xq(i)+u)/dsqrt(gam) - sum=sum+wq(i)*x**l*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) - enddo - int_prod_bessel_large=sum*factor -end - -!! Calculation of -!! -!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) -!! -!! M_n(x) modified spherical bessel function -!! -double precision function int_prod_bessel_loc(l,gam,n,a) - implicit none - integer n,k,l,kcp - double precision gam,a - double precision int,intold,coef_nk,crochet,dble_fact, fact, pi, expo - double precision :: f_0, f_k - logical done - - pi=dacos(-1.d0) - intold=-1.d0 - done=.false. - int=0 - - ! Int f_0 - coef_nk=1.d0/dble_fact( n+n+1 ) - expo=0.5d0*dfloat(n+l+1) - crochet=dble_fact(n+l-1)/(gam+gam)**expo - if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(0.5d0*pi) - - f_0 = coef_nk*a**n*crochet - - k=0 - - f_k=f_0 - do while (.not.done) - - int=int+f_k - -! f_k = f_k*(a**2*(2*(k+1)+n+l-1)) / (2*(k+1)*(2*(n+k+1)+1)*2*gam) - f_k = f_k*(a*a*dble(k+k+1+n+l)) / (dble((k+k+2)*(4*(n+k+1)+2))*gam) - - if(dabs(int-intold).lt.1d-15)then - done=.true. - else - k=k+1 - intold=int - endif - - enddo - - int_prod_bessel_loc=int -end - - double precision function int_prod_bessel_num(l,gam,n,m,a,b) - implicit none - integer n,m,l,i,npoints - double precision gam,a,b - double precision sum,dx,x,bessel_mod - sum=0.d0 - npoints=20000 - dx=30.d0/npoints - do i=1,npoints - x=(i-1)*dx+0.5d0*dx - sum=sum+dx*x**l*dexp(-gam*x**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m) - enddo - int_prod_bessel_num=sum - end - - -! l,m : Y(l,m) parameters -! c(3) : pseudopotential center -! a(3) : Atomic Orbital center -! n_a(3) : Powers of x,y,z in the Atomic Orbital -! g_a : Atomic Orbital exponent -! r : Distance between the Atomic Orbital center and the considered point -double precision function ylm_orb(l,m,c,a,n_a,g_a,r) -implicit none -integer lmax_max -integer l,m -double precision a(3),g_a,c(3) -double precision prod,binom_func,accu,bigI,ylm,bessel_mod -double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal -integer ntotA,mu,k1,k2,k3,lambda -integer n_a(3) -double precision y, f1, f2 -double precision, allocatable :: array_coefs_A(:,:) - -ac2=(a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2 -ac=dsqrt(ac2) -arg=g_a*(ac2+r*r) -fourpi=4.d0*dacos(-1.d0) -factor=fourpi*dexp(-arg) -areal=2.d0*g_a*ac -ntotA=n_a(1)+n_a(2)+n_a(3) - - -if(ac.eq.0.d0)then - ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) - return -else - - theta_AC0=dacos( (a(3)-c(3))/ac ) - phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) - - allocate (array_coefs_A(0:ntotA,3)) - do k1=0,n_a(1) - array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)*r**(k1) - enddo - do k2=0,n_a(2) - array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)*r**(k2) - enddo - do k3=0,n_a(3) - array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)*r**(k3) - enddo - - accu=0.d0 - do lambda=0,l+ntotA - do mu=-lambda,lambda - y = ylm(lambda,mu,theta_AC0,phi_AC0) - if (y == 0.d0) then - cycle - endif - do k3=0,n_a(3) - f1 = y*array_coefs_A(k3,3) - if (f1 == 0.d0) cycle - do k2=0,n_a(2) - f2 = f1*array_coefs_A(k2,2) - if (f2 == 0.d0) cycle - do k1=0,n_a(1) - prod=f2*array_coefs_A(k1,1)*bigI(lambda,mu,l,m,k1,k2,k3) - if (prod == 0.d0) then - cycle - endif - if (areal*r < 100.d0) then ! overflow! - accu=accu+prod*bessel_mod(areal*r,lambda) - endif - enddo - enddo - enddo - enddo - enddo - ylm_orb=factor*accu - deallocate (array_coefs_A) - return -endif - -end diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index 339156db..89315a92 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -3,7 +3,7 @@ logical function ao_one_e_integral_zero(i,k) integer, intent(in) :: i,k ao_one_e_integral_zero = .False. - if (.not.((io_ao_integrals_overlap/='None').or.is_periodic.or.use_cgtos)) then + if((io_ao_integrals_overlap=='None').and.(.not.use_cgtos)) then if (ao_overlap_abs(i,k) < ao_one_e_integrals_threshold) then ao_one_e_integral_zero = .True. return diff --git a/src/ao_one_e_ints/spread_dipole_ao.irp.f b/src/ao_one_e_ints/spread_dipole_ao.irp.f index 86469a3f..d00ea871 100644 --- a/src/ao_one_e_ints/spread_dipole_ao.irp.f +++ b/src/ao_one_e_ints/spread_dipole_ao.irp.f @@ -1,376 +1,68 @@ - BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)] + implicit none BEGIN_DOC ! * array of the integrals of AO_i * x^2 AO_j - ! + END_DOC + call ao_cart_to_ao_basis(ao_cart_spread_x, ao_cart_num, ao_spread_x, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)] + implicit none + BEGIN_DOC ! * array of the integrals of AO_i * y^2 AO_j - ! + END_DOC + call ao_cart_to_ao_basis(ao_cart_spread_y, ao_cart_num, ao_spread_y, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)] + implicit none + BEGIN_DOC ! * array of the integrals of AO_i * z^2 AO_j END_DOC - implicit none - integer :: i,j,n,l - double precision :: f, tmp - integer :: dim1 - double precision :: overlap, overlap_x, overlap_y, overlap_z - double precision :: alpha, beta - double precision :: A_center(3), B_center(3) - integer :: power_A(3), power_B(3) - double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z - dim1=500 - lower_exp_val = 40.d0 - ao_spread_x= 0.d0 - ao_spread_y= 0.d0 - ao_spread_z= 0.d0 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_x,overlap_y, overlap_z, overlap, & - !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - accu_x = 0.d0 - accu_y = 0.d0 - accu_z = 0.d0 - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i) - beta = ao_expo_ordered_transp(l,i) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) - accu_x += c*tmp*overlap_y*overlap_z - call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) - accu_y += c*tmp*overlap_x*overlap_z - call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) - accu_z += c*tmp*overlap_y*overlap_x - enddo - enddo - ao_spread_x(i,j) = accu_x - ao_spread_y(i,j) = accu_y - ao_spread_z(i,j) = accu_z - enddo - enddo - !$OMP END PARALLEL DO + call ao_cart_to_ao_basis(ao_cart_spread_z, ao_cart_num, ao_spread_z, ao_num) END_PROVIDER - - BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)] BEGIN_DOC ! * array of the integrals of AO_i * x AO_j - ! + END_DOC + implicit none + call ao_cart_to_ao_basis(ao_cart_dipole_x, ao_cart_num, ao_dipole_x, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)] + BEGIN_DOC ! * array of the integrals of AO_i * y AO_j - ! + END_DOC + implicit none + call ao_cart_to_ao_basis(ao_cart_dipole_y, ao_cart_num, ao_dipole_y, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)] + BEGIN_DOC ! * array of the integrals of AO_i * z AO_j END_DOC implicit none - integer :: i,j,n,l - double precision :: f, tmp - integer :: dim1 - double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z - double precision :: alpha, beta - double precision :: A_center(3), B_center(3) - integer :: power_A(3), power_B(3) - double precision :: lower_exp_val, dx, c - dim1=500 - lower_exp_val = 40.d0 - ao_dipole_x= 0.d0 - ao_dipole_y= 0.d0 - ao_dipole_z= 0.d0 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_x,overlap_y, overlap_z, overlap, & - !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - accu_x = 0.d0 - accu_y = 0.d0 - accu_z = 0.d0 - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - - call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) - accu_x = accu_x + c*tmp*overlap_y*overlap_z - call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) - accu_y = accu_y + c*tmp*overlap_x*overlap_z - call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) - accu_z = accu_z + c*tmp*overlap_y*overlap_x - enddo - enddo - ao_dipole_x(i,j) = accu_x - ao_dipole_y(i,j) = accu_y - ao_dipole_z(i,j) = accu_z - enddo - enddo - !$OMP END PARALLEL DO + call ao_cart_to_ao_basis(ao_cart_dipole_z, ao_cart_num, ao_dipole_z, ao_num) END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)] BEGIN_DOC ! * array of the integrals of AO_i * d/dx AO_j - ! + END_DOC + implicit none + call ao_cart_to_ao_basis(ao_cart_deriv_1_x, ao_cart_num, ao_deriv_1_x, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)] + BEGIN_DOC ! * array of the integrals of AO_i * d/dy AO_j - ! + END_DOC + implicit none + call ao_cart_to_ao_basis(ao_cart_deriv_1_y, ao_cart_num, ao_deriv_1_y, ao_num) + END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)] + BEGIN_DOC ! * array of the integrals of AO_i * d/dz AO_j END_DOC implicit none - integer :: i,j,n,l - double precision :: f, tmp - integer :: dim1 - double precision :: overlap, overlap_x, overlap_y, overlap_z - double precision :: alpha, beta - double precision :: A_center(3), B_center(3) - integer :: power_A(3), power_B(3) - double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z - integer :: i_component - dim1=500 - lower_exp_val = 40.d0 - ao_deriv_1_x= 0.d0 - ao_deriv_1_y= 0.d0 - ao_deriv_1_z= 0.d0 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(A_center,B_center,power_A,power_B,& - !$OMP overlap_x,overlap_y, overlap_z, overlap, & - !$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) & - !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & - !$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & - !$OMP ao_expo_ordered_transp,dim1,lower_exp_val) - do j=1,ao_num - A_center(1) = nucl_coord( ao_nucl(j), 1 ) - A_center(2) = nucl_coord( ao_nucl(j), 2 ) - A_center(3) = nucl_coord( ao_nucl(j), 3 ) - power_A(1) = ao_power( j, 1 ) - power_A(2) = ao_power( j, 2 ) - power_A(3) = ao_power( j, 3 ) - do i= 1,ao_num - B_center(1) = nucl_coord( ao_nucl(i), 1 ) - B_center(2) = nucl_coord( ao_nucl(i), 2 ) - B_center(3) = nucl_coord( ao_nucl(i), 3 ) - power_B(1) = ao_power( i, 1 ) - power_B(2) = ao_power( i, 2 ) - power_B(3) = ao_power( i, 3 ) - accu_x = 0.d0 - accu_y = 0.d0 - accu_z = 0.d0 - do n = 1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(n,j) - do l = 1, ao_prim_num(i) - beta = ao_expo_ordered_transp(l,i) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - c = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(n,j) - i_component = 1 - call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) - accu_x += c*(tmp*overlap_y*overlap_z) - i_component = 2 - call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) - accu_y += c*(tmp*overlap_x*overlap_z) - i_component = 3 - call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) - accu_z += c*(tmp*overlap_y*overlap_x) - enddo - enddo - ao_deriv_1_x(i,j) = accu_x - ao_deriv_1_y(i,j) = accu_y - ao_deriv_1_z(i,j) = accu_z - enddo - enddo - !$OMP END PARALLEL DO + call ao_cart_to_ao_basis(ao_cart_deriv_1_z, ao_cart_num, ao_deriv_1_z, ao_num) END_PROVIDER - - - - subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) - BEGIN_DOC -! Computes the following integral : -! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x^2 ] -! needed for the dipole and those things - END_DOC - implicit none - integer :: i,j,k,l - integer,intent(in) :: power_A,power_B - double precision, intent(in) :: lower_exp_val - double precision,intent(in) :: A_center, B_center,alpha,beta - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho - double precision :: P_center,pouet_timy - if(power_A.lt.0.or.power_B.lt.0)then - overlap_x = 0.d0 - dx = 0.d0 - return - endif - p = alpha + beta - p_inv= 1.d0/p - rho = alpha * beta * p_inv - dist = (A_center - B_center)*(A_center - B_center) - P_center = (alpha * A_center + beta * B_center) * p_inv - factor = dexp(-rho * dist) - if(factor.lt.0.000001d0)then -! print*,'factor = ',factor - dx = 0.d0 - overlap_x = 0.d0 - return - endif - pouet_timy = dsqrt(lower_exp_val/p) - x_min = P_center - pouet_timy - x_max = P_center + pouet_timy - domain = x_max-x_min - dx = domain/dble(nx) - overlap_x = 0.d0 - x = x_min - do i = 1, nx - x += dx - overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x * x - enddo - overlap_x *= factor * dx - - end - - - subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) -! compute the following integral : -! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] -! needed for the dipole and those things - implicit none - integer :: i,j,k,l - integer,intent(in) :: power_A,power_B - double precision, intent(in) :: lower_exp_val - double precision,intent(in) :: A_center, B_center,alpha,beta - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho - double precision :: P_center - if(power_A.lt.0.or.power_B.lt.0)then - overlap_x = 0.d0 - dx = 0.d0 - return - endif - p = alpha + beta - p_inv= 1.d0/p - rho = alpha * beta * p_inv - dist = (A_center - B_center)*(A_center - B_center) - P_center = (alpha * A_center + beta * B_center) * p_inv - factor = dexp(-rho * dist) - if(power_B == 0 .and. power_A ==0)then - double precision :: F_integral - overlap_x = P_center * F_integral(0,p) * factor - dx = 0.d0 - return - endif - double precision :: pouet_timy - - pouet_timy = dsqrt(lower_exp_val/p) - x_min = P_center - pouet_timy - x_max = P_center + pouet_timy - domain = x_max-x_min - dx = domain/dble(nx) - overlap_x = 0.d0 - x = x_min - do i = 1, nx - x += dx - overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x - enddo - overlap_x *= factor * dx - - end - - subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx) - implicit none - integer :: i,j,k,l - integer,intent(in) :: power_A(3),power_B(3),i_component - double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: overlap_first, overlap_second -! computes : = (a_x_i - 2 alpha ) - - call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx) - call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx) - overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second) - end - - subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) - implicit none -! compute the following integral : -! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] - integer :: i,j,k,l - integer,intent(in) :: power_A,power_B - double precision, intent(in) :: lower_exp_val - double precision,intent(in) :: A_center, B_center,alpha,beta - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho - double precision :: P_center,pouet_timy - if(power_A.lt.0.or.power_B.lt.0)then - overlap_x = 0.d0 - dx = 0.d0 - return - endif - p = alpha + beta - p_inv= 1.d0/p - rho = alpha * beta * p_inv - dist = (A_center - B_center)*(A_center - B_center) - P_center = (alpha * A_center + beta * B_center) * p_inv - factor = dexp(-rho * dist) - if(factor.lt.0.000001d0)then - dx = 0.d0 - overlap_x = 0.d0 - return - endif - - pouet_timy = dsqrt(lower_exp_val/p) - x_min = P_center - pouet_timy - x_max = P_center + pouet_timy - domain = x_max-x_min - dx = domain/dble(nx) - overlap_x = 0.d0 - x = x_min - do i = 1, nx - x += dx - overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) - enddo - overlap_x *= factor * dx - end - diff --git a/src/ao_one_e_ints/test_ao_one_ints.irp.f b/src/ao_one_e_ints/test_ao_one_ints.irp.f new file mode 100644 index 00000000..8c39b352 --- /dev/null +++ b/src/ao_one_e_ints/test_ao_one_ints.irp.f @@ -0,0 +1,27 @@ + +BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] + implicit none + BEGIN_DOC +! overlap matrix of the ao_ortho_canonical. +! Expected to be the Identity + END_DOC + integer :: i,j,k,l + double precision :: c + do j=1, ao_ortho_canonical_num + do i=1, ao_ortho_canonical_num + ao_ortho_canonical_overlap(i,j) = 0.d0 + enddo + enddo + do j=1, ao_ortho_canonical_num + do k=1, ao_num + c = 0.d0 + do l=1, ao_num + c += ao_ortho_canonical_coef(l,j) * ao_overlap(l,k) + enddo + do i=1, ao_ortho_canonical_num + ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c + enddo + enddo + enddo +END_PROVIDER +