diff --git a/config/ifort.cfg b/config/ifort.cfg index 2b2fe0a2..100b87af 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -31,7 +31,7 @@ OPENMP : 1 ; Append OpenMP flags # -ftz : Flushes denormal results to zero # [OPT] -FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g +FCFLAGS : -axSSE4.2,AVX -O2 -ip -ftz -g # Profiling flags ################# diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 27a79b4f..e16909bb 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -24,9 +24,9 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num_align,ao_prim_num BEGIN_DOC ! Coefficients including the AO normalization END_DOC - double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3) + double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3), c integer :: l, powA(3), nz - integer :: i,j + integer :: i,j,k nz=100 C_A(1) = 0.d0 C_A(2) = 0.d0 @@ -39,6 +39,17 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num_align,ao_prim_num call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) enddo + ! Normalization of the contracted basis functions + norm = 0.d0 + do j=1,ao_prim_num(i) + do k=1,ao_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) + enddo + enddo + do j=1,ao_prim_num(i) + ao_coef_normalized(i,j) = ao_coef_normalized(i,j)/sqrt(norm) + enddo enddo END_PROVIDER diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index b1f1dba0..dd19f9d4 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -40,24 +40,22 @@ double precision function ao_bielec_integral(i,j,k,l) L_center(p) = nucl_coord(num_l,p) enddo + double precision :: coef1, coef2, coef3, coef4 + double precision :: p_inv,q_inv + double precision :: general_primitive_integral + do p = 1, ao_prim_num(i) - double precision :: coef1 coef1 = ao_coef_normalized_ordered_transp(p,i) do q = 1, ao_prim_num(j) - double precision :: coef2 coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - double precision :: p_inv,q_inv call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & I_power,J_power,I_center,J_center,dim1) p_inv = 1.d0/pp do r = 1, ao_prim_num(k) - double precision :: coef3 coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) do s = 1, ao_prim_num(l) - double precision :: coef4 coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - double precision :: general_primitive_integral call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & K_power,L_power,K_center,L_center,dim1) @@ -65,7 +63,7 @@ double precision function ao_bielec_integral(i,j,k,l) integral = general_primitive_integral(dim1, & P_new,P_center,fact_p,pp,p_inv,iorder_p, & Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_bielec_integral += coef4 * integral + ao_bielec_integral = ao_bielec_integral + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -94,7 +92,7 @@ double precision function ao_bielec_integral(i,j,k,l) I_power(1),J_power(1),K_power(1),L_power(1), & I_power(2),J_power(2),K_power(2),L_power(2), & I_power(3),J_power(3),K_power(3),L_power(3)) - ao_bielec_integral += coef4 * integral + ao_bielec_integral = ao_bielec_integral + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -479,11 +477,11 @@ double precision function general_primitive_integral(dim, & enddo n_Ix = 0 do ix = 0, iorder_p(1) - if (abs(P_new(ix,1)) < 1.d-8) cycle + if (abs(P_new(ix,1)) < ao_integrals_threshold) cycle a = P_new(ix,1) do jx = 0, iorder_q(1) d = a*Q_new(jx,1) - if (abs(d) < 1.d-8) cycle + if (abs(d) < ao_integrals_threshold) cycle !DEC$ FORCEINLINE call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx) !DEC$ FORCEINLINE @@ -500,11 +498,11 @@ double precision function general_primitive_integral(dim, & enddo n_Iy = 0 do iy = 0, iorder_p(2) - if (abs(P_new(iy,2)) > 1.d-8) then + if (abs(P_new(iy,2)) > ao_integrals_threshold) then b = P_new(iy,2) do jy = 0, iorder_q(2) e = b*Q_new(jy,2) - if (abs(e) < 1.d-8) cycle + if (abs(e) < ao_integrals_threshold) cycle !DEC$ FORCEINLINE call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny) !DEC$ FORCEINLINE @@ -522,11 +520,11 @@ double precision function general_primitive_integral(dim, & enddo n_Iz = 0 do iz = 0, iorder_p(3) - if (abs(P_new(iz,3)) > 1.d-8) then + if (abs(P_new(iz,3)) > ao_integrals_threshold) then c = P_new(iz,3) do jz = 0, iorder_q(3) f = c*Q_new(jz,3) - if (abs(f) < 1.d-8) cycle + if (abs(f) < ao_integrals_threshold) cycle !DEC$ FORCEINLINE call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz) !DEC$ FORCEINLINE diff --git a/src/Integrals_Monoelec/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f index eadc0b72..3c23b458 100644 --- a/src/Integrals_Monoelec/pot_ao_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_ints.irp.f @@ -171,7 +171,7 @@ include 'Utils/constants.include.F' enddo const_factor = dist*rho const = p * dist_integral - if(const_factor.ge.80.d0)then + if(const_factor > 80.d0)then NAI_pol_mult = 0.d0 return endif @@ -375,10 +375,10 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) Y(ix) = 0.d0 enddo call I_x2_pol_mult_mono_elec(c-1,R1x,R1xp,R2x,X,nx,n_pt_in) - do ix=0,nx - X(ix) *= c - enddo - call multiply_poly(X,nx,R2x,2,d,nd) + do ix=0,nx + X(ix) *= dble(c) + enddo + call multiply_poly(X,nx,R2x,2,d,nd) ny=0 call I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,Y,ny,n_pt_in) call multiply_poly(Y,ny,R1x,2,d,nd) @@ -390,10 +390,10 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) nx = 0 call I_x1_pol_mult_mono_elec(a-2,c,R1x,R1xp,R2x,X,nx,n_pt_in) ! print*,'nx a-2,c= ',nx - do ix=0,nx - X(ix) *= a-1 - enddo - call multiply_poly(X,nx,R2x,2,d,nd) + do ix=0,nx + X(ix) *= dble(a-1) + enddo + call multiply_poly(X,nx,R2x,2,d,nd) ! print*,'nd out = ',nd nx = nd @@ -403,7 +403,7 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) call I_x1_pol_mult_mono_elec(a-1,c-1,R1x,R1xp,R2x,X,nx,n_pt_in) ! print*,'nx a-1,c-1 = ',nx do ix=0,nx - X(ix) *= c + X(ix) *= dble(c) enddo call multiply_poly(X,nx,R2x,2,d,nd) ny=0 @@ -444,7 +444,7 @@ recursive subroutine I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,d,nd,dim) call I_x1_pol_mult_mono_elec(0,c-2,R1x,R1xp,R2x,X,nx,dim) ! print*,'nx 0,c-2 = ',nx do ix=0,nx - X(ix) *= c-1 + X(ix) *= dble(c-1) enddo call multiply_poly(X,nx,R2x,2,d,nd) ! print*,'nd = ',nd diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 59693cdb..f3efaa4c 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -68,7 +68,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu c = 0.d0 if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& - < 1.d-10) then + < ao_integrals_threshold) then cycle endif do k = 1, nucl_num @@ -165,10 +165,10 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu c = 0.d0 if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& - < 1.d-10) then + < ao_integrals_threshold) then cycle endif - + do k = 1, nucl_num double precision :: Z Z = nucl_charge(k) diff --git a/src/Integrals_Monoelec/spread_dipole_ao.irp.f b/src/Integrals_Monoelec/spread_dipole_ao.irp.f index d7aa738a..5611ec7f 100644 --- a/src/Integrals_Monoelec/spread_dipole_ao.irp.f +++ b/src/Integrals_Monoelec/spread_dipole_ao.irp.f @@ -55,11 +55,11 @@ 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) + 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) + 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) + accu_z += c*tmp*overlap_y*overlap_x enddo enddo ao_spread_x(i,j) = accu_x @@ -130,11 +130,11 @@ 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) + 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) + 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) + accu_z = accu_z + c*tmp*overlap_y*overlap_x enddo enddo ao_dipole_x(i,j) = accu_x diff --git a/src/Integrals_Monoelec/var_pt2_ratio.irp.f b/src/Integrals_Monoelec/var_pt2_ratio.irp.f deleted file mode 100644 index e69de29b..00000000 diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 746b5f13..bcfd43b8 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -83,13 +83,22 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) D(i) = 1.d0/dsqrt(D(i)) else m = i-1 + print *, 'Removed Linear dependencies below:', 1.d0/D(m) exit endif enddo do i=m+1,n + print *, D(i) D(i) = 0.d0 enddo + do i=1,m + if ( D(i) >= 1.d5 ) then + print *, 'Warning: Basis set may have linear dependence problems' + endif + enddo + + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S_half,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j) diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F index 632dc50b..b55880c3 100644 --- a/src/Utils/constants.include.F +++ b/src/Utils/constants.include.F @@ -1,4 +1,4 @@ -integer, parameter :: max_dim = 255 +integer, parameter :: max_dim = 511 integer, parameter :: SIMD_vector = 32 double precision, parameter :: pi = dacos(-1.d0) diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f index e9c3f9ab..ade130f7 100644 --- a/src/Utils/integration.irp.f +++ b/src/Utils/integration.irp.f @@ -78,7 +78,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, !DEC$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) - if (fact_k < 1.d-8) then + if (fact_k < ao_integrals_threshold) then fact_k = 0.d0 return endif @@ -210,7 +210,7 @@ subroutine gaussian_product(a,xa,b,xb,k,p,xp) xab(3) = xa(3)-xb(3) ab = ab*p_inv k = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3)) - if (k > 20.d0) then + if (k > 40.d0) then k=0.d0 return endif @@ -249,7 +249,7 @@ subroutine gaussian_product_x(a,xa,b,xb,k,p,xp) xab = xa-xb ab = ab*p_inv k = ab*xab*xab - if (k > 20.d0) then + if (k > 40.d0) then k=0.d0 return endif @@ -580,7 +580,7 @@ double precision function rint_large_n(n,rho) enddo t2=0.d0 do l=0,k - t2=t2+(-1.d0)**l/fact(l+1)/fact(k-l) + t2=t2+(-1.d0)**l/(fact(l+1)*fact(k-l)) enddo alpha_k=t2*fact(k+1)*fact(k)*(-1.d0)**k alpha_k= alpha_k/t1 diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index fd15054e..1b74f430 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -150,19 +150,19 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& integer :: i do i = 1,iorder_p(1) - overlap_x += P_new(i,1) * F_integral_tab(i) + overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) enddo call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1)) overlap_x *= fact_p do i = 1,iorder_p(2) - overlap_y += P_new(i,2) * F_integral_tab(i) + overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) enddo call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2)) overlap_y *= fact_p do i = 1,iorder_p(3) - overlap_z += P_new(i,3) * F_integral_tab(i) + overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) enddo call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3)) overlap_z *= fact_p