From 8703a06c4faf7616e5ea276f27f1e1da84bf0f1d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 17 Mar 2017 19:05:30 +0100 Subject: [PATCH] Fixed pseudos --- src/AO_Basis/aos.irp.f | 5 +- .../pot_ao_pseudo_ints.irp.f | 68 ++++++++++--------- src/Integrals_Monoelec/pseudopot.f90 | 46 +++++-------- 3 files changed, 55 insertions(+), 64 deletions(-) diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 0938d3bd..f0f03fab 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -182,7 +182,7 @@ integer function ao_power_index(nx,ny,nz) end -BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)] +BEGIN_PROVIDER [ character*(128), l_to_charater, (0:7)] BEGIN_DOC ! character corresponding to the "L" value of an AO orbital END_DOC @@ -192,6 +192,9 @@ BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)] l_to_charater(2)='D' l_to_charater(3)='F' l_to_charater(4)='G' + l_to_charater(5)='H' + l_to_charater(6)='I' + l_to_charater(7)='J' END_PROVIDER diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index fa855113..7cfd6f9f 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -51,22 +51,23 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu print*, 'Providing the nuclear electron pseudo integrals (local)' call wall_time(wall_1) + wall_0 = wall_1 call cpu_time(cpu_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,n_pt_in, & -! !$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_integral_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() -! !$OMP DO SCHEDULE (guided) + !$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 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_integral_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() + !$OMP DO SCHEDULE (guided) do j = 1, ao_num @@ -120,8 +121,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu endif enddo -! !$OMP END DO -! !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER @@ -148,23 +149,24 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu print*, 'Providing the nuclear electron pseudo integrals (non-local)' call wall_time(wall_1) + wall_0 = wall_1 call cpu_time(cpu_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,n_pt_in, & -! !$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_integral_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() -! -! !$OMP DO SCHEDULE (guided) + !$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 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_integral_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() + + !$OMP DO SCHEDULE (guided) +! do j = 1, ao_num num_A = ao_nucl(j) @@ -217,12 +219,12 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu endif endif enddo -! -! !$OMP END DO -! -! !$OMP END PARALLEL - - + + !$OMP END DO + + !$OMP END PARALLEL + + END_PROVIDER BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ] diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index d77b3ca0..a69aa42d 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -15,14 +15,10 @@ double precision function Vps & implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3) -integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) -integer lmax,kmax,n_kl(kmax_max,0:lmax_max) -double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) -integer klocmax_max -parameter (klocmax_max=10) -integer klocmax,n_k(klocmax_max) -double precision v_k(klocmax_max),dz_k(klocmax_max) +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) & @@ -36,13 +32,10 @@ double precision function Vps_num & implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3),rmax -integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) -integer lmax,kmax,n_kl(kmax_max,0:lmax_max) -double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) -integer klocmax_max;parameter (klocmax_max=10) -integer klocmax,n_k(klocmax_max) -double precision v_k(klocmax_max),dz_k(klocmax_max) +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 @@ -54,11 +47,9 @@ 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_max -parameter (klocmax_max=10) integer klocmax -double precision v_k(klocmax_max),dz_k(klocmax_max) -integer n_k(klocmax_max) +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) @@ -705,12 +696,9 @@ end 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_max,lmax_max,ntot_max -parameter (klocmax_max=10,lmax_max=2) -parameter (ntot_max=10) integer klocmax -double precision v_k(klocmax_max),dz_k(klocmax_max),crochet,bigA -integer n_k(klocmax_max) +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 @@ -719,6 +707,7 @@ 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) @@ -755,8 +744,8 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg dreal=2.d0*d2 - allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max)) - allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) + 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 @@ -2111,9 +2100,7 @@ end ! 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,ntot_max -parameter (lmax_max=2) -parameter (ntot_max=14) +integer lmax_max integer l,m double precision a(3),g_a,c(3) double precision prod,binom_func,accu,bigI,ylm,bessel_mod @@ -2131,7 +2118,6 @@ factor=fourpi*dexp(-arg) areal=2.d0*g_a*ac ntotA=n_a(1)+n_a(2)+n_a(3) -if(ntotA.gt.ntot_max)stop 'increase ntot_max' 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))