From 48a001a353ab8f3bbb20545c06f4c67f39d7ed74 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 Jan 2012 09:37:21 +0100 Subject: [PATCH] Bug --- src/ao.irp.f | 35 ++++++++++++++++++++++++++--------- src/ao_oneD_point.irp.f | 14 +++++++------- src/eplf_function.irp.f | 28 ++++++++++++++-------------- 3 files changed, 47 insertions(+), 30 deletions(-) diff --git a/src/ao.irp.f b/src/ao.irp.f index 8bf774e..d9466c7 100644 --- a/src/ao.irp.f +++ b/src/ao.irp.f @@ -118,20 +118,20 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ] END_PROVIDER - BEGIN_PROVIDER [ real, ao_expo, (ao_num,ao_prim_num_max) ] -&BEGIN_PROVIDER [ real, ao_coef, (ao_num,ao_prim_num_max) ] + BEGIN_PROVIDER [ real, ao_expo_transp, (ao_num,ao_prim_num_max) ] +&BEGIN_PROVIDER [ real, ao_coef_transp, (ao_num,ao_prim_num_max) ] implicit none BEGIN_DOC ! Exponents and coefficients of the atomic orbitals END_DOC - ao_expo = 0. - call get_ao_basis_ao_expo(ao_expo) + ao_expo_transp = 0. + call get_ao_basis_ao_expo(ao_expo_transp) integer :: i,j do i=1,ao_num 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 write(message,'(A,I6,A,I6,A)') 'Exponent ',j,' of contracted gaussian ',i,' is < 0' call abrt(irp_here,message) @@ -139,8 +139,8 @@ END_PROVIDER enddo enddo - ao_coef = 0. - call get_ao_basis_ao_coef(ao_coef) + ao_coef_transp = 0. + call get_ao_basis_ao_coef(ao_coef_transp) ! Normalization of the AO coefficients ! ------------------------------------ @@ -152,8 +152,25 @@ END_PROVIDER pow(1) = ao_power(i,1) pow(2) = ao_power(i,2) pow(3) = ao_power(i,3) - norm = goverlap(ao_expo(i,j),ao_expo(i,j),pow) - ao_coef(i,j) = ao_coef(i,j)/sqrt(norm) + norm = goverlap(ao_expo_transp(i,j),ao_expo_transp(i,j),pow) + 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 diff --git a/src/ao_oneD_point.irp.f b/src/ao_oneD_point.irp.f index 18965c0..3da66f3 100644 --- a/src/ao_oneD_point.irp.f +++ b/src/ao_oneD_point.irp.f @@ -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) do k=1,ao_prim_num_max 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 @@ -46,7 +46,7 @@ BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ] enddo do k=1,ao_prim_num_max 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 @@ -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 i=1,ao_num 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 @@ -85,7 +85,7 @@ BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ] enddo do k=1,ao_prim_num_max 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 @@ -101,8 +101,8 @@ BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_num,ao_prim_num_max) ] integer :: i, k, l do k=1,ao_prim_num_max do i=1,ao_num - ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo(i,k) * & - ( 4.*ao_expo(i,k)*point_nucl_dist_2(ao_nucl(i)) - 6. ) + ao_oneD_prim_lapl_p(i,k) = ao_oneD_prim_p(i,k) * ao_expo_transp(i,k) * & + ( 4.*ao_expo_transp(i,k)*point_nucl_dist_2(ao_nucl(i)) - 6. ) enddo enddo @@ -122,7 +122,7 @@ BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ] enddo do k=1,ao_prim_num_max 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 diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 99399e7..ea244e3 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -198,31 +198,31 @@ END_PROVIDER ! ao_eplf_integral_numeric = 0.d0 ! do q=1,ao_prim_num(j) ! 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 = & ! ao_eplf_integral_primitive_oneD_numeric( & -! ao_expo(i,p), & +! ao_expo(p,i), & ! nucl_coord(ao_nucl(i),1), & ! ao_power(i,1), & -! ao_expo(j,q), & +! ao_expo(q,j), & ! nucl_coord(ao_nucl(j),1), & ! ao_power(j,1), & ! gmma, & ! center(1)) * & ! ao_eplf_integral_primitive_oneD_numeric( & -! ao_expo(i,p), & +! ao_expo(p,i), & ! nucl_coord(ao_nucl(i),2), & ! ao_power(i,2), & -! ao_expo(j,q), & +! ao_expo(q,j), & ! nucl_coord(ao_nucl(j),2), & ! ao_power(j,2), & ! gmma, & ! center(2)) * & ! ao_eplf_integral_primitive_oneD_numeric( & -! ao_expo(i,p), & +! ao_expo(p,i), & ! nucl_coord(ao_nucl(i),3), & ! ao_power(i,3), & -! ao_expo(j,q), & +! ao_expo(q,j), & ! nucl_coord(ao_nucl(j),3), & ! ao_power(j,3), & ! gmma, & @@ -314,10 +314,10 @@ double precision function ao_eplf_integral(i,j,gmma,center) do p=1,ao_prim_num(i) integral = & ao_eplf_integral_primitive_oneD( & - ao_expo(i,p), & + ao_expo(p,i), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & - ao_expo(j,q), & + ao_expo(q,j), & nucl_coord(ao_nucl(j),1), & ao_power(j,1), & gmma, & @@ -325,10 +325,10 @@ double precision function ao_eplf_integral(i,j,gmma,center) if (integral /= 0.d0) then integral *= & ao_eplf_integral_primitive_oneD( & - ao_expo(i,p), & + ao_expo(p,i), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & - ao_expo(j,q), & + ao_expo(q,j), & nucl_coord(ao_nucl(j),2), & ao_power(j,2), & gmma, & @@ -336,15 +336,15 @@ double precision function ao_eplf_integral(i,j,gmma,center) if (integral /= 0.d0) then integral *= & ao_eplf_integral_primitive_oneD( & - ao_expo(i,p), & + ao_expo(p,i), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & - ao_expo(j,q), & + ao_expo(q,j), & nucl_coord(ao_nucl(j),3), & ao_power(j,3), & gmma, & 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 enddo