LCPQ
/
eplf
mirror of https://gitlab.com/scemama/eplf
10
0
Fork 0
Code Issues Releases Wiki Activity

Bug

This commit is contained in:
Anthony Scemama 2012-01-26 09:37:21 +01:00
parent 57f3d60e52
commit 48a001a353
3 changed files with 47 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
src/ao.irp.f
View File

@ -118,20 +118,20 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ real, ao_expo, (ao_num,ao_prim_num_max) ] BEGIN_PROVIDER [ real, ao_expo_transp, (ao_num,ao_prim_num_max) ]
&BEGIN_PROVIDER [ real, ao_coef, (ao_num,ao_prim_num_max) ] &BEGIN_PROVIDER [ real, ao_coef_transp, (ao_num,ao_prim_num_max) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Exponents and coefficients of the atomic orbitals ! Exponents and coefficients of the atomic orbitals
END_DOC END_DOC
ao_expo = 0. ao_expo_transp = 0.
call get_ao_basis_ao_expo(ao_expo) call get_ao_basis_ao_expo(ao_expo_transp)
integer :: i,j integer :: i,j
do i=1,ao_num do i=1,ao_num
do j=1,ao_prim_num(i) do j=1,ao_prim_num(i)
if (ao_expo(i,j) <= 0.) then if (ao_expo_transp(i,j) <= 0.) then
character*(80) :: message character*(80) :: message
write(message,'(A,I6,A,I6,A)') 'Exponent ',j,' of contracted gaussian ',i,' is < 0' write(message,'(A,I6,A,I6,A)') 'Exponent ',j,' of contracted gaussian ',i,' is < 0'
call abrt(irp_here,message) call abrt(irp_here,message)
@ -139,8 +139,8 @@ END_PROVIDER
enddo enddo
enddo enddo
ao_coef = 0. ao_coef_transp = 0.
call get_ao_basis_ao_coef(ao_coef) call get_ao_basis_ao_coef(ao_coef_transp)
! Normalization of the AO coefficients ! Normalization of the AO coefficients
! ------------------------------------ ! ------------------------------------
@ -152,8 +152,25 @@ END_PROVIDER
pow(1) = ao_power(i,1) pow(1) = ao_power(i,1)
pow(2) = ao_power(i,2) pow(2) = ao_power(i,2)
pow(3) = ao_power(i,3) pow(3) = ao_power(i,3)
norm = goverlap(ao_expo(i,j),ao_expo(i,j),pow) norm = goverlap(ao_expo_transp(i,j),ao_expo_transp(i,j),pow)
ao_coef(i,j) = ao_coef(i,j)/sqrt(norm) ao_coef_transp(i,j) = ao_coef_transp(i,j)/sqrt(norm)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, ao_expo, (ao_prim_num_max,ao_num) ]
&BEGIN_PROVIDER [ real, ao_coef, (ao_prim_num_max,ao_num) ]
implicit none
BEGIN_DOC
! Exponents and coefficients of the atomic orbitals
END_DOC
integer :: i,j
do i=1,ao_num
do j=1,ao_prim_num(i)
ao_coef(j,i) = ao_coef_transp(i,j)
ao_expo(j,i) = ao_expo_transp(i,j)
enddo enddo
enddo enddo

14
src/ao_oneD_point.irp.f
View File

@ -17,7 +17,7 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_num,ao_prim_num_max) ]
! Compute exp(-alpha*r) or exp(-alpha*r^2) ! Compute exp(-alpha*r) or exp(-alpha*r^2)
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
ao_oneD_prim_p(i,k) = exp(-ao_oneD_prim_p(i,k)*ao_expo(i,k)) ao_oneD_prim_p(i,k) = exp(-ao_oneD_prim_p(i,k)*ao_expo_transp(i,k))
enddo enddo
enddo enddo
@ -46,7 +46,7 @@ BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ]
enddo enddo
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(i,k)*ao_oneD_prim_p(i,k) ao_oneD_p(i) = ao_oneD_p(i) + ao_coef_transp(i,k)*ao_oneD_prim_p(i,k)
enddo enddo
enddo enddo
@ -65,7 +65,7 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_num,ao_prim_num_max,3) ]
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
factor = -2.*point_nucl_dist_vec(ao_nucl(i),l) factor = -2.*point_nucl_dist_vec(ao_nucl(i),l)
ao_oneD_prim_grad_p(i,k,l) = factor*ao_expo(i,k)*ao_oneD_prim_p(i,k) ao_oneD_prim_grad_p(i,k,l) = factor*ao_expo_transp(i,k)*ao_oneD_prim_p(i,k)
enddo enddo
enddo enddo
enddo enddo
@ -85,7 +85,7 @@ BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ]
enddo enddo
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(i,k)*ao_oneD_prim_grad_p(i,k,l) ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef_transp(i,k)*ao_oneD_prim_grad_p(i,k,l)
enddo enddo
enddo enddo
enddo enddo
@ -101,8 +101,8 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_num,ao_prim_num_max) ]
integer :: i, k, l integer :: i, k, l
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo(i,k) * & ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo_transp(i,k) * &
( 4.*ao_expo(i,k)*point_nucl_dist_2(ao_nucl(i)) - 6. ) ( 4.*ao_expo_transp(i,k)*point_nucl_dist_2(ao_nucl(i)) - 6. )
enddo enddo
enddo enddo
@ -122,7 +122,7 @@ BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ]
enddo enddo
do k=1,ao_prim_num_max do k=1,ao_prim_num_max
do i=1,ao_num do i=1,ao_num
ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(i,k)*ao_oneD_prim_lapl_p(i,k) ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef_transp(i,k)*ao_oneD_prim_lapl_p(i,k)
enddo enddo
enddo enddo

28
src/eplf_function.irp.f
View File

@ -198,31 +198,31 @@ END_PROVIDER
! ao_eplf_integral_numeric = 0.d0 ! ao_eplf_integral_numeric = 0.d0
! do q=1,ao_prim_num(j) ! do q=1,ao_prim_num(j)
! do p=1,ao_prim_num(i) ! do p=1,ao_prim_num(i)
! c = ao_coef(i,p)*ao_coef(j,q) ! c = ao_coef(p,i)*ao_coef(q,j)
! integral = & ! integral = &
! ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
! ao_expo(i,p), & ! ao_expo(p,i), &
! nucl_coord(ao_nucl(i),1), & ! nucl_coord(ao_nucl(i),1), &
! ao_power(i,1), & ! ao_power(i,1), &
! ao_expo(j,q), & ! ao_expo(q,j), &
! nucl_coord(ao_nucl(j),1), & ! nucl_coord(ao_nucl(j),1), &
! ao_power(j,1), & ! ao_power(j,1), &
! gmma, & ! gmma, &
! center(1)) * & ! center(1)) * &
! ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
! ao_expo(i,p), & ! ao_expo(p,i), &
! nucl_coord(ao_nucl(i),2), & ! nucl_coord(ao_nucl(i),2), &
! ao_power(i,2), & ! ao_power(i,2), &
! ao_expo(j,q), & ! ao_expo(q,j), &
! nucl_coord(ao_nucl(j),2), & ! nucl_coord(ao_nucl(j),2), &
! ao_power(j,2), & ! ao_power(j,2), &
! gmma, & ! gmma, &
! center(2)) * & ! center(2)) * &
! ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
! ao_expo(i,p), & ! ao_expo(p,i), &
! nucl_coord(ao_nucl(i),3), & ! nucl_coord(ao_nucl(i),3), &
! ao_power(i,3), & ! ao_power(i,3), &
! ao_expo(j,q), & ! ao_expo(q,j), &
! nucl_coord(ao_nucl(j),3), & ! nucl_coord(ao_nucl(j),3), &
! ao_power(j,3), & ! ao_power(j,3), &
! gmma, & ! gmma, &
@ -314,10 +314,10 @@ double precision function ao_eplf_integral(i,j,gmma,center)
do p=1,ao_prim_num(i) do p=1,ao_prim_num(i)
integral = & integral = &
ao_eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(i,p), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),1), & nucl_coord(ao_nucl(i),1), &
ao_power(i,1), & ao_power(i,1), &
ao_expo(j,q), & ao_expo(q,j), &
nucl_coord(ao_nucl(j),1), & nucl_coord(ao_nucl(j),1), &
ao_power(j,1), & ao_power(j,1), &
gmma, & gmma, &
@ -325,10 +325,10 @@ double precision function ao_eplf_integral(i,j,gmma,center)
if (integral /= 0.d0) then if (integral /= 0.d0) then
integral *= & integral *= &
ao_eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(i,p), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),2), & nucl_coord(ao_nucl(i),2), &
ao_power(i,2), & ao_power(i,2), &
ao_expo(j,q), & ao_expo(q,j), &
nucl_coord(ao_nucl(j),2), & nucl_coord(ao_nucl(j),2), &
ao_power(j,2), & ao_power(j,2), &
gmma, & gmma, &
@ -336,15 +336,15 @@ double precision function ao_eplf_integral(i,j,gmma,center)
if (integral /= 0.d0) then if (integral /= 0.d0) then
integral *= & integral *= &
ao_eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(i,p), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),3), & nucl_coord(ao_nucl(i),3), &
ao_power(i,3), & ao_power(i,3), &
ao_expo(j,q), & ao_expo(q,j), &
nucl_coord(ao_nucl(j),3), & nucl_coord(ao_nucl(j),3), &
ao_power(j,3), & ao_power(j,3), &
gmma, & gmma, &
center(3)) center(3))
buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q) buffer(p) += integral*ao_coef(p,i)*ao_coef(q,j)
endif endif
endif endif
enddo enddo