LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00
Code Issues Releases Wiki Activity

Using normalized atomic basis functions now

Browse Source
This commit is contained in:
Anthony Scemama 2015-12-11 14:32:41 +01:00
parent 204984bdd8
commit 13857ccebe
11 changed files with 63 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -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
#################

15
src/AO_Basis/aos.irp.f
View File

@ -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

26
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -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

22
src/Integrals_Monoelec/pot_ao_ints.irp.f
View File

@ -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

6
src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
View File

@ -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)

12
src/Integrals_Monoelec/spread_dipole_ao.irp.f
View File

@ -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

0
src/Integrals_Monoelec/var_pt2_ratio.irp.f
View File

9
src/Utils/LinearAlgebra.irp.f
View File

@ -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)

2
src/Utils/constants.include.F
View File

@ -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)

8
src/Utils/integration.irp.f
View File

@ -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

6
src/Utils/one_e_integration.irp.f
View File

@ -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