diff --git a/src/FCIdump/README.rst b/src/FCIdump/README.rst index 8a467b16..faf69203 100644 --- a/src/FCIdump/README.rst +++ b/src/FCIdump/README.rst @@ -10,9 +10,6 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fcidump `_ - Undocumented - Needed Modules diff --git a/src/Pseudo_integrals/README.rst b/src/Pseudo_integrals/README.rst index 08076556..4f98b220 100644 --- a/src/Pseudo_integrals/README.rst +++ b/src/Pseudo_integrals/README.rst @@ -9,6 +9,12 @@ Documentation .. NEEDED_MODULES file. `ao_nucl_elec_integral_pseudo `_ + Undocumented + +`ao_nucl_elec_integral_pseudo_local `_ + Local component of the pseudopotential + +`ao_nucl_elec_integral_pseudo_non_local `_ interaction nuclear electron diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90 index 4eb9d2ae..0cea5f66 100644 --- a/src/Pseudo_integrals/int.f90 +++ b/src/Pseudo_integrals/int.f90 @@ -1003,11 +1003,11 @@ end double precision function crochet(n,g) implicit none integer n -double precision g,pi,dble_fact,expo -pi=dacos(-1.d0) +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)/(2.d0*g)**expo -if(mod(n,2).eq.0)crochet=crochet*dsqrt(pi/2.d0) +if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2 end !! @@ -2043,10 +2043,10 @@ double precision function int_prod_bessel_loc(l,gam,n,a) int=0 ! Int f_0 - coef_nk=1.d0/dble_fact(2*n+1) + coef_nk=1.d0/dble_fact( n+n+1 ) expo=0.5d0*dfloat(n+l+1) crochet=dble_fact(n+l-1)/(2.d0*gam)**expo - if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(pi/2.d0) + if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(0.5d0*pi) f_0 = coef_nk*a**n*crochet diff --git a/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f b/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f index 37dc7591..ca39b402 100644 --- a/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f +++ b/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f @@ -1,8 +1,17 @@ - BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num)] - BEGIN_DOC -! interaction nuclear electron - END_DOC +BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num) ] implicit none +! Nucleus-pseudopotential interaction + END_DOC + write(output_monoints,*) 'Providing the pseudopotential integrals ' + ao_nucl_elec_integral_pseudo = ao_nucl_elec_integral_pseudo_local + ao_nucl_elec_integral_pseudo_non_local +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo_local, (ao_num_align,ao_num)] + implicit none + BEGIN_DOC +! Local component of the pseudopotential + END_DOC double precision :: alpha, beta, gama, delta integer :: num_A,num_B double precision :: A_center(3),B_center(3),C_center(3) @@ -13,13 +22,8 @@ double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 integer :: thread_num - ao_nucl_elec_integral_pseudo = 0.d0 - - ! - ! | _ _ _. | - ! |_ (_) (_ (_| | - ! - !! Parameters of the local part of pseudo: + PROVIDE output_monoints + ao_nucl_elec_integral_pseudo_local = 0.d0 integer klocmax integer, allocatable :: n_k(:,:) @@ -40,41 +44,6 @@ allocate(n_k_dump(1:klocmax), v_k_dump(1:klocmax), dz_k_dump(1:klocmax)) - ! - ! |\ | _ ._ | _ _ _. | - ! | \| (_) | | | (_) (_ (_| | - ! - !! Parameters of non local part of pseudo: - - integer :: kmax,lmax - integer, allocatable :: n_kl(:,:,:) - double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) - - call ezfio_get_pseudo_integrals_lmaxpo(lmax) - call ezfio_get_pseudo_integrals_kmax(kmax) - !lmax plus one -> lmax - lmax = lmax - 1 - - allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) - - call ezfio_get_pseudo_integrals_n_kl(n_kl) - call ezfio_get_pseudo_integrals_v_kl(v_kl) - call ezfio_get_pseudo_integrals_dz_kl(dz_kl) - - - !! Dump array - integer, allocatable :: n_kl_dump(:,:) - double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) - - allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) - - ! _ - ! / _. | _ | - ! \_ (_| | (_ |_| | - ! - - write(output_monoints,*) 'Providing the nuclear electron pseudo integrals ' - call wall_time(wall_1) call cpu_time(cpu_1) @@ -82,11 +51,11 @@ !$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,n_pt_in, & - !$OMP v_k_dump,n_k_dump, dz_k_dump, n_kl_dump, v_kl_dump, dz_kl_dump, & + !$OMP v_k_dump,n_k_dump, dz_k_dump, & !$OMP wall_0,wall_2,thread_num, output_monoints) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & - !$OMP ao_nucl_elec_integral_pseudo,nucl_num,nucl_charge, & - !$OMP klocmax,lmax,kmax,v_k,n_k, dz_k, n_kl, v_kl, dz_kl, & + !$OMP ao_nucl_elec_integral_pseudo_local,nucl_num,nucl_charge, & + !$OMP klocmax,v_k,n_k, dz_k, & !$OMP wall_1) !$OMP DO SCHEDULE (guided) @@ -124,15 +93,8 @@ c = c + Vloc(klocmax, v_k_dump,n_k_dump, dz_k_dump, & A_center,power_A,alpha,B_center,power_B,beta,C_center) - - n_kl_dump = n_kl(k,1:kmax,0:lmax) - v_kl_dump = v_kl(k,1:kmax,0:lmax) - dz_kl_dump = dz_kl(k,1:kmax,0:lmax) - - c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) - enddo - ao_nucl_elec_integral_pseudo(i,j) = ao_nucl_elec_integral_pseudo(i,j) + & + ao_nucl_elec_integral_pseudo_local(i,j) = ao_nucl_elec_integral_pseudo_local(i,j) + & ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c enddo enddo @@ -152,16 +114,134 @@ !$OMP END PARALLEL -! _ -! | \ _ _. | | _ _ _. _|_ _ -! |_/ (/_ (_| | | (_) (_ (_| |_ (/_ -! - deallocate(n_k,v_k, dz_k) deallocate(n_k_dump,v_k_dump, dz_k_dump) + + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo_non_local, (ao_num_align,ao_num)] + BEGIN_DOC +! interaction nuclear electron + END_DOC + implicit none + 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,n_pt_in,m + double precision :: Vloc, Vpseudo + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + integer :: thread_num + + PROVIDE output_monoints + ao_nucl_elec_integral_pseudo_non_local = 0.d0 + + integer :: kmax,lmax + integer, allocatable :: n_kl(:,:,:) + double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) + + call ezfio_get_pseudo_integrals_lmaxpo(lmax) + call ezfio_get_pseudo_integrals_kmax(kmax) + !lmax plus one -> lmax + lmax = lmax - 1 + + allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) + + call ezfio_get_pseudo_integrals_n_kl(n_kl) + call ezfio_get_pseudo_integrals_v_kl(v_kl) + call ezfio_get_pseudo_integrals_dz_kl(dz_kl) + + + !! Dump array + integer, allocatable :: n_kl_dump(:,:) + double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) + + allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) + + ! _ + ! / _. | _ | + ! \_ (_| | (_ |_| | + ! + + call wall_time(wall_1) + call cpu_time(cpu_1) + + !$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,n_pt_in, & + !$OMP n_kl_dump, v_kl_dump, dz_kl_dump, & + !$OMP wall_0,wall_2,thread_num, output_monoints) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & + !$OMP ao_nucl_elec_integral_pseudo_non_local,nucl_num,nucl_charge, & + !$OMP lmax,kmax,v_kl, dz_kl, n_kl, & + !$OMP 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) + beta = ao_expo_ordered_transp(m,i) + double precision :: c + c = 0.d0 + + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + n_kl_dump = n_kl(k,1:kmax,0:lmax) + v_kl_dump = v_kl(k,1:kmax,0:lmax) + dz_kl_dump = dz_kl(k,1:kmax,0:lmax) + + c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + ao_nucl_elec_integral_pseudo_non_local(i,j) = ao_nucl_elec_integral_pseudo_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 + write(output_monoints,*) 100.*float(j)/float(ao_num), '% in ', & + wall_2-wall_1, 's' + endif + endif + enddo + + !$OMP END DO + !$OMP END PARALLEL + + deallocate(n_kl,v_kl, dz_kl) deallocate(n_kl_dump,v_kl_dump, dz_kl_dump) END_PROVIDER + + + + diff --git a/src/Utils/README.rst b/src/Utils/README.rst index ff16f3ed..b36fb9c6 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -242,7 +242,7 @@ Documentation `align_double `_ Compute 1st dimension such that it is aligned for vectorization. -`approx_dble `_ +`approx_dble `_ Undocumented `binom `_ @@ -257,10 +257,16 @@ Documentation `binom_transp `_ Binomial coefficients -`dble_fact `_ +`dble_fact `_ + Undocumented + +`dble_fact_odd `_ n!! -`dble_logfact `_ +`dble_fact_peer `_ + n!! + +`dble_logfact `_ n!! `fact `_ @@ -269,29 +275,29 @@ Documentation `fact_inv `_ 1/n! -`inv_int `_ +`inv_int `_ 1/i `logfact `_ n! -`normalize `_ +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads -`u_dot_u `_ +`u_dot_u `_ Compute -`u_dot_v `_ +`u_dot_v `_ Compute -`wall_time `_ +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 7a78ad59..705f4f7b 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -138,21 +138,21 @@ END_PROVIDER double precision function dble_fact(n) implicit none integer :: n - double precision :: dble_fact_peer, dble_fact_odd + double precision :: dble_fact_even, dble_fact_odd dble_fact = 1.d0 if(n.lt.0) return if(iand(n,1).eq.0)then - dble_fact = dble_fact_peer(n) + dble_fact = dble_fact_even(n) else dble_fact= dble_fact_odd(n) endif end function -double precision function dble_fact_peer(n) result(fact2) +double precision function dble_fact_even(n) result(fact2) implicit none BEGIN_DOC ! n!!