From 2a6497119578c5d815940f9455fac8cbfcda5e9b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 14 May 2009 17:48:27 +0200 Subject: [PATCH] d_up_up is negative for open shell => bug --- ao_axis_point.irp.f | 107 ++++++++++++++ ao_full_point.irp.f | 51 +++++++ ao_oneD_point.irp.f | 128 +++++++++++++++++ electrons.irp.f | 54 +++++++ eplf.irp.f | 337 ++++++++++++++++++++++++++++++++++++++++++++ mo.irp.f | 136 ++++++++++++++++++ point.irp.f | 48 +++++++ test_1d.irp.f | 18 +++ 8 files changed, 879 insertions(+) create mode 100644 ao_axis_point.irp.f create mode 100644 ao_full_point.irp.f create mode 100644 ao_oneD_point.irp.f create mode 100644 electrons.irp.f create mode 100644 eplf.irp.f create mode 100644 mo.irp.f create mode 100644 point.irp.f create mode 100644 test_1d.irp.f diff --git a/ao_axis_point.irp.f b/ao_axis_point.irp.f new file mode 100644 index 0000000..2e93c19 --- /dev/null +++ b/ao_axis_point.irp.f @@ -0,0 +1,107 @@ +BEGIN_PROVIDER [ real, ao_axis_power_p, (-2:ao_power_max,3,nucl_num) ] + implicit none + + BEGIN_DOC +! Evaluation of power of x, y, z at the current point for each +! nucleus. Negative power -> 0. + END_DOC + + integer :: i,k,l + do i=1,nucl_num + do l=1,3 + ao_axis_power_p(-2,l,i) = 0. + ao_axis_power_p(-1,l,i) = 0. + ao_axis_power_p(0,l,i) = 0. + ao_axis_power_p(0,l,i) = 1. + do k=1,ao_power_max_nucl(i,l) + ao_axis_power_p(k,l,i) = point_nucl_dist_vec(i,l)*ao_axis_power_p(k-1,l,i) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_axis_p, (ao_num) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Cartesian polynomial part of the atomic orbitals. + END_DOC + integer :: i + + do i=1,ao_num + ao_axis_p(i) & + = ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_axis_grad_p, (ao_num,3) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Gradients of the cartesian polynomial part of the atomic orbitals. + END_DOC + integer :: i, l + real:: real_of_int(-1:10) + data real_of_int /0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./ + + do i=1,ao_num + ao_axis_grad_p(i,1) = real_of_int(ao_power(i,1)) & + * ao_axis_power_p( ao_power(i,1)-1, 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) + enddo + + do i=1,ao_num + ao_axis_grad_p(i,2) = real_of_int(ao_power(i,2)) & + * ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2)-1, 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) + enddo + + do i=1,ao_num + ao_axis_grad_p(i,3) = real_of_int(ao_power(i,3)) & + * ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3)-1, 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_axis_lapl_p, (ao_num) ] + implicit none + include 'types.F' + BEGIN_DOC +! Laplacian of the cartesian atomic orbitals + END_DOC + integer :: i, j, l + + do i=1,ao_num + real:: real_of_int(-2:10) + data real_of_int /0.,0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./ + + ao_axis_lapl_p(i) & + = real_of_int(ao_power(i,1)) & + * real_of_int(ao_power(i,1)-1) & + * ao_axis_power_p( ao_power(i,1)-2, 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) & + + real_of_int(ao_power(i,2)) & + * real_of_int(ao_power(i,2)-1) & + * ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2)-2, 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) & + + real_of_int(ao_power(i,3)) & + * real_of_int(ao_power(i,3)-1) & + * ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) & + * ao_axis_power_p( ao_power(i,3)-2, 3 , ao_nucl(i) ) + enddo + +END_PROVIDER + + diff --git a/ao_full_point.irp.f b/ao_full_point.irp.f new file mode 100644 index 0000000..f7f74d7 --- /dev/null +++ b/ao_full_point.irp.f @@ -0,0 +1,51 @@ +BEGIN_PROVIDER [ real, ao_value_p, (ao_num) ] + implicit none + BEGIN_DOC +! Values of the atomic orbitals + END_DOC + integer :: i + do i=1,ao_num + ao_value_p(i) = ao_oneD_p(i) * ao_axis_p(i) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_grad_p, (ao_num,3) ] + implicit none + include 'types.F' + BEGIN_DOC + ! Gradients of the atomic orbitals + END_DOC + + integer :: i,l + do l=1,3 + do i=1,ao_num + ao_grad_p(i,l) = ao_oneD_p(i) * ao_axis_grad_p(i,l) + enddo + do i=1,ao_num + ao_grad_p(i,l) = ao_grad_p(i,l) + ao_oneD_grad_p(i,l) * ao_axis_p(i) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_lapl_p, (ao_num) ] + implicit none + include 'types.F' + BEGIN_DOC +! Laplacian of the atomic orbitals + END_DOC + integer :: i,l + do i=1,ao_num + ao_lapl_p(i) = ao_oneD_p(i) * ao_axis_lapl_p(i) + enddo + do i=1,ao_num + ao_lapl_p(i) = ao_lapl_p(i) + ao_oneD_lapl_p(i) * ao_axis_p(i) + enddo + do l=1,3 + do i=1,ao_num + ao_lapl_p(i) = ao_lapl_p(i) + 2.*ao_oneD_grad_p(i,l) * ao_axis_grad_p(i,l) + enddo + enddo + +END_PROVIDER + diff --git a/ao_oneD_point.irp.f b/ao_oneD_point.irp.f new file mode 100644 index 0000000..5897810 --- /dev/null +++ b/ao_oneD_point.irp.f @@ -0,0 +1,128 @@ +BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_prim_num_max,ao_num) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Exponentials of the primitive AOs + END_DOC + integer :: i, k + real:: r2, rtemp + +! Compute alpha*r or alpha*r^2 + do i=1,ao_num + r2 = point_nucl_dist_2(ao_nucl(i)) + do k=1,ao_prim_num(i) + ao_oneD_prim_p(k,i) = r2 + enddo + enddo + + ! Compute exp(-alpha*r) or exp(-alpha*r^2) + do i=1,ao_num + do k=1,ao_prim_num(i) + ao_oneD_prim_p(k,i) = exp(-ao_oneD_prim_p(k,i)*ao_expo(k,i)) + enddo + ! Cut below 1.d-12 + do k=1,ao_prim_num(i) + if ( abs(ao_oneD_prim_p(k,i)) < 1.e-12 ) then + ao_oneD_prim_p(k,i) = 0. + endif + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ] + implicit none + include 'types.F' + + BEGIN_DOC +! One-dimensional component of the AOs + END_DOC + + integer :: i, k + + do i=1,ao_num + ao_oneD_p(i) = 0. + do k=1,ao_prim_num(i) + ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(k,i)*ao_oneD_prim_p(k,i) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_prim_num_max,ao_num,3) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Gradients of the one-dimensional component of the primitive AOs + END_DOC + integer :: i, k, l + real:: factor + do l=1,3 + do i=1,ao_num + factor = -2.*point_nucl_dist_vec(ao_nucl(i),l) + do k=1,ao_prim_num(i) + ao_oneD_prim_grad_p(k,i,l) = factor*ao_expo(k,i)*ao_oneD_prim_p(k,i) + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Gradients of the one-dimensional component of the AOs + END_DOC + integer :: i, k, l + do l=1,3 + do i=1,ao_num + ao_oneD_grad_p(i,l) = 0. + do k=1,ao_prim_num(i) + ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(k,i)*ao_oneD_prim_grad_p(k,i,l) + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_prim_num_max,ao_num) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Laplacian of the one-dimensional component of the primitive AOs + END_DOC + integer :: i, k, l + do i=1,ao_num + do k=1,ao_prim_num(i) + ao_oneD_prim_lapl_p(k,i) = ao_oneD_prim_p(k,i) * ao_expo(k,i) * & + ( 4.*ao_expo(k,i)*point_nucl_dist_2(ao_nucl(i)) - 6. ) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ] + implicit none + include 'types.F' + + BEGIN_DOC +! Laplacian of the one-dimensional component of the AOs + END_DOC + + integer :: i, k + + do i=1,ao_num + ao_oneD_lapl_p(i) = 0. + do k=1,ao_prim_num(i) + ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(k,i)*ao_oneD_prim_lapl_p(k,i) + enddo + enddo + +END_PROVIDER + diff --git a/electrons.irp.f b/electrons.irp.f new file mode 100644 index 0000000..651831e --- /dev/null +++ b/electrons.irp.f @@ -0,0 +1,54 @@ +BEGIN_PROVIDER [ integer, elec_alpha_num ] + + BEGIN_DOC +! Number of alpha electrons + END_DOC + + implicit none + elec_alpha_num = -1 +!$OMP CRITICAL (qcio_critical) + call qcio_get_system_num_alpha(elec_alpha_num) +!$OMP END CRITICAL (qcio_critical) + ASSERT (elec_alpha_num > 0) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, elec_beta_num ] + + BEGIN_DOC +! Number of beta electrons + END_DOC + + implicit none + elec_beta_num = -1 +!$OMP CRITICAL (qcio_critical) + call qcio_get_system_num_beta(elec_beta_num) +!$OMP END CRITICAL (qcio_critical) + ASSERT (elec_beta_num >= 0) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, elec_num ] + + BEGIN_DOC +! Number of electrons + END_DOC + implicit none + + elec_num = elec_alpha_num + elec_beta_num + + ASSERT ( elec_num > 0 ) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, elec_num_2, (2) ] + + BEGIN_DOC +! Number of alpha and beta electrons in an array + END_DOC + + elec_num_2(1) = elec_alpha_num + elec_num_2(2) = elec_beta_num + +END_PROVIDER + diff --git a/eplf.irp.f b/eplf.irp.f new file mode 100644 index 0000000..596738c --- /dev/null +++ b/eplf.irp.f @@ -0,0 +1,337 @@ +BEGIN_PROVIDER [ real, eplf_gamma ] + implicit none + BEGIN_DOC +! Value of the gaussian for the EPLF + END_DOC + eplf_gamma = 10. +END_PROVIDER + +BEGIN_PROVIDER [ double precision, eplf_integral_matrix, (ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Array of all the for EPLF + END_DOC + integer :: i, j + double precision :: eplf_integral + do i=1,ao_num + do j=i,ao_num + eplf_integral_matrix(j,i) = eplf_integral(j,i,eplf_gamma,point) + eplf_integral_matrix(i,j) = eplf_integral_matrix(j,i) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, eplf_up_up ] +&BEGIN_PROVIDER [ double precision, eplf_up_dn ] + implicit none + BEGIN_DOC +! Value of the d_upup and d_updn quantities needed for EPLF + END_DOC + + integer :: i, j, k, l, m, n + + eplf_up_up = 0. + eplf_up_dn = 0. + + PROVIDE mo_coef_transp + + do m=1,ao_num + do n=1,ao_num + double precision :: ao_mn + ao_mn = eplf_integral_matrix(m,n) + + if (ao_mn /= 0.d0) then + + do k=1,ao_num + do l=1,ao_num + double precision :: ao_klmn + ao_klmn = ao_mn*ao_value_p(k)*ao_value_p(l) + + if (ao_klmn /= 0.d0) then + + do j=1,elec_beta_num + if (mo_coef_transp(j,n) /= 0.d0) then + do i=1,elec_beta_num + eplf_up_up = eplf_up_up + 2.d0*ao_klmn* & + mo_coef_transp(i,k)*mo_coef_transp(j,n) * & + (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & + mo_coef_transp(i,m)*mo_coef_transp(j,l) ) + enddo + endif + enddo + + do j=1+elec_beta_num, elec_alpha_num + if (mo_coef_transp(j,n) /= 0.d0) then + do i=1,elec_beta_num + eplf_up_up = eplf_up_up + 2.d0*ao_klmn* & + mo_coef_transp(i,k)*mo_coef_transp(j,n) * & + (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & + mo_coef_transp(i,m)*mo_coef_transp(j,l) ) + enddo + do i=1+elec_beta_num,elec_alpha_num + eplf_up_up = eplf_up_up + ao_klmn* & + mo_coef_transp(i,k)*mo_coef_transp(j,n) * & + (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & + mo_coef_transp(i,m)*mo_coef_transp(j,l) ) + enddo + endif + enddo + + do j=1,elec_beta_num + if ( (mo_coef_transp(j,n) /= 0.d0).and. & + (mo_coef_transp(j,m) /= 0.d0) ) then + do i=1,elec_alpha_num + eplf_up_dn = eplf_up_dn + 2.d0*ao_klmn* & + mo_coef_transp(i,k)*mo_coef_transp(j,n) * & + mo_coef_transp(i,l)*mo_coef_transp(j,m) + enddo + endif + enddo + + endif + + enddo + enddo + + endif + + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ real, eplf_value ] + implicit none + BEGIN_DOC +! Value of the EPLF at the current point. + END_DOC + double precision :: aa, ab + + aa = eplf_up_up + ab = eplf_up_dn + aa = min(1.d0,aa) + ab = min(1.d0,ab) + aa = max(1.d-30,aa) + ab = max(1.d-30,ab) + aa = -dlog(aa)/eplf_gamma + ab = -dlog(ab)/eplf_gamma + aa = dsqrt(aa) + ab = dsqrt(ab) + + eplf_value = (aa-ab)/(aa+ab) + +END_PROVIDER + + +double precision function eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr) + implicit none + include 'constants.F' + + real, intent(in) :: a,b,gmma ! Exponents + real, intent(in) :: xa,xb,xr ! Centers + integer, intent(in) :: i,j ! Powers of xa and xb + integer,parameter :: Npoints=1000 + real :: x, xmin, xmax, dx + + ASSERT (a>0.) + ASSERT (b>0.) + ASSERT (i>=0) + ASSERT (j>=0) + + xmin = min(xa,xb) + xmax = max(xa,xb) + xmin = min(xmin,xr) - 10. + xmax = max(xmax,xr) + 10. + dx = (xmax-xmin)/real(Npoints) + + real :: dtemp + dtemp = 0. + x = xmin + integer :: k + do k=1,Npoints + dtemp = dtemp + & + (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) + x = x+dx + enddo + eplf_integral_primitive_oneD_numeric = dtemp*dx + +end function + +double precision function eplf_integral_numeric(i,j,gmma,center) + implicit none + integer, intent(in) :: i, j + integer :: p,q,k + double precision :: integral(ao_prim_num_max,ao_prim_num_max) + double precision :: eplf_integral_primitive_oneD_numeric + real :: gmma, center(3) + + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + integral(p,q) = & + eplf_integral_primitive_oneD_numeric( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),1), & + ao_power(i,1), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),1), & + ao_power(j,1), & + gmma, & + center(1)) * & + eplf_integral_primitive_oneD_numeric( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),2), & + ao_power(i,2), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),2), & + ao_power(j,2), & + gmma, & + center(2)) * & + eplf_integral_primitive_oneD_numeric( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),3), & + ao_power(i,3), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),3), & + ao_power(j,3), & + gmma, & + center(3)) + enddo + enddo + + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j) + enddo + enddo + + eplf_integral_numeric = 0. + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + eplf_integral_numeric = eplf_integral_numeric + integral(p,q) + enddo + enddo + +end function + +double precision function eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) + implicit none + include 'constants.F' + + real, intent(in) :: a,b,gmma ! Exponents + real, intent(in) :: xa,xb,xr ! Centers + integer, intent(in) :: i,j ! Powers of xa and xb + integer :: ii, jj, kk, ll + real :: xp1,xp + real :: p1,p + double precision :: S(0:i+1,0:j) + double precision :: inv_p, di(max(i,j)), dj(j), c, c1 + + ASSERT (a>0.) + ASSERT (b>0.) + ASSERT (i>=0) + ASSERT (j>=0) + + ! Gaussian product + call gaussian_product(a,xa,b,xb,c1,p1,xp1) + call gaussian_product(p1,xp1,gmma,xr,c,p,xp) + inv_p = 1.d0/p + S(0,0) = dsqrt(pi*inv_p)*c*c1 + + ! Obara-Saika recursion + + if (i>0) then + S(1,0) = (xp-xa) * S(0,0) + endif + if (j>0) then + S(0,1) = (xp-xb) * S(0,0) + endif + + do ii=1,max(i,j) + di(ii) = 0.5d0*inv_p*dble(ii) + enddo + + if (i>1) then + do ii=1,i-1 + S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0) + enddo + endif + + if (j>1) then + do jj=1,j-1 + S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1) + enddo + endif + + do jj=1,j + S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1) + do ii=2,i + S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1) + enddo + enddo + + eplf_integral_primitive_oneD = S(i,j) + +end function + +double precision function eplf_integral(i,j,gmma,center) + implicit none + integer, intent(in) :: i, j + integer :: p,q,k + double precision :: integral(ao_prim_num_max,ao_prim_num_max) + double precision :: eplf_integral_primitive_oneD + real :: gmma, center(3) + + ASSERT(i>0) + ASSERT(j>0) + ASSERT(i<=ao_num) + ASSERT(j<=ao_num) + + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + integral(p,q) = & + eplf_integral_primitive_oneD( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),1), & + ao_power(i,1), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),1), & + ao_power(j,1), & + gmma, & + center(1)) * & + eplf_integral_primitive_oneD( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),2), & + ao_power(i,2), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),2), & + ao_power(j,2), & + gmma, & + center(2)) * & + eplf_integral_primitive_oneD( & + ao_expo(p,i), & + nucl_coord(ao_nucl(i),3), & + ao_power(i,3), & + ao_expo(q,j), & + nucl_coord(ao_nucl(j),3), & + ao_power(j,3), & + gmma, & + center(3)) + enddo + enddo + + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j) + enddo + enddo + + eplf_integral = 0. + do q=1,ao_prim_num(j) + do p=1,ao_prim_num(i) + eplf_integral = eplf_integral + integral(p,q) + enddo + enddo + +end function + diff --git a/mo.irp.f b/mo.irp.f new file mode 100644 index 0000000..72dbc3a --- /dev/null +++ b/mo.irp.f @@ -0,0 +1,136 @@ +BEGIN_PROVIDER [ integer, mo_closed_num ] + implicit none + + BEGIN_DOC +! Number of closed shell Molecular orbitals + END_DOC + +!$OMP CRITICAL (qcio_critical) + call qcio_get_mo_num_closed(mo_closed_num) +!$OMP END CRITICAL (qcio_critical) + ASSERT (mo_closed_num >= 0) + ASSERT (mo_closed_num <= elec_alpha_num) + ASSERT (mo_closed_num <= elec_beta_num) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, mo_active_num ] + implicit none + + BEGIN_DOC +! Number of active Molecular orbitals + END_DOC + +!$OMP CRITICAL (qcio_critical) + call qcio_get_mo_num_active(mo_active_num) +!$OMP END CRITICAL (qcio_critical) + ASSERT (mo_active_num >= 0) + +END_PROVIDER + +BEGIN_PROVIDER [ integer, mo_num ] + implicit none + + BEGIN_DOC +! Number of Molecular orbitals + END_DOC + + mo_num = mo_closed_num + mo_active_num + ASSERT(mo_num > 0) + +END_PROVIDER + + +BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ] + implicit none + + BEGIN_DOC +! Molecular orbital coefficients + END_DOC + + integer :: i, j + double precision :: buffer(ao_num,mo_tot_num) + +!$OMP CRITICAL (qcio_critical) + call qcio_get_mo_matrix(buffer) +!$OMP END CRITICAL (qcio_critical) + do j=1,mo_num + do i=1,ao_num + mo_coef(i,j) = buffer(i,j) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ logical, mo_non_zero, (mo_num,ao_num) ] + implicit none + + BEGIN_DOC + ! Values where the mo coefficients are /= 0. + END_DOC + + integer :: i, j, k + + do j=1,ao_num + do k=1,mo_num + mo_non_zero(k,j) = mo_coef_transp(j,k) /= 0. + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num,ao_num) ] + implicit none + + BEGIN_DOC + ! Transposed array of mo coefficients + END_DOC + + integer :: i, j, k + + do j=1,ao_num + do k=1,mo_num + mo_coef_transp(k,j) = mo_coef(j,k) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ logical, mo_is_closed, (mo_num) ] +&BEGIN_PROVIDER [ logical, mo_is_active, (mo_num) ] + implicit none + + BEGIN_DOC +! mo_is_closed : True if mo(i) is a closed-shell +! mo_is_active : True if mo(i) is an active orbital + END_DOC + + character :: buffer(mo_tot_num) + + integer :: i + do i=1,mo_num + if ( buffer(i) == 'c' ) then + mo_is_closed(i) = .True. + mo_is_active(i) = .False. + else if ( buffer(i) == 'a' ) then + mo_is_closed(i) = .False. + mo_is_active(i) = .True. + else + mo_is_closed(i) = .False. + mo_is_active(i) = .False. + endif + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, mo_tot_num ] + + BEGIN_DOC +! Total number of MOs in the QCIO file + END_DOC +!$OMP CRITICAL (qcio_critical) + call qcio_get_mo_num_orb_tot(mo_tot_num) +!$OMP END CRITICAL (qcio_critical) + ASSERT (mo_tot_num > 0) + +END_PROVIDER + diff --git a/point.irp.f b/point.irp.f new file mode 100644 index 0000000..a368d3f --- /dev/null +++ b/point.irp.f @@ -0,0 +1,48 @@ +BEGIN_PROVIDER [ real, point, (3) ] + implicit none + BEGIN_DOC +! Coordinates of the current point + END_DOC + point(1) = 0. + point(2) = 0. + point(3) = 0. +END_PROVIDER + +BEGIN_PROVIDER [ real, point_nucl_dist_vec, (nucl_num,3) ] + implicit none + BEGIN_DOC +! Distance vector between the current point and the nuclei + END_DOC + integer :: k + do k=1,nucl_num + point_nucl_dist_vec(k,1) = point(1)-nucl_coord(k,1) + point_nucl_dist_vec(k,2) = point(2)-nucl_coord(k,2) + point_nucl_dist_vec(k,3) = point(3)-nucl_coord(k,3) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ real, point_nucl_dist_2, (nucl_num) ] + implicit none + BEGIN_DOC +! Square of the distance between the current point and the nuclei + END_DOC + integer :: k,l + do k=1,nucl_num + point_nucl_dist_2(k) = 0. + do l=1,3 + point_nucl_dist_2(k) = point_nucl_dist_2(k) + point_nucl_dist_vec(k,l)**2 + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ real, point_nucl_dist, (nucl_num) ] + implicit none + BEGIN_DOC +! Distance between the current point and the nuclei + END_DOC + integer :: k + do k=1,nucl_num + point_nucl_dist(k) = sqrt(point_nucl_dist_2(k)) + enddo +END_PROVIDER + diff --git a/test_1d.irp.f b/test_1d.irp.f new file mode 100644 index 0000000..660c188 --- /dev/null +++ b/test_1d.irp.f @@ -0,0 +1,18 @@ +program debug + PROVIDE ao_prim_num_max + read(*,*) eplf_gamma + TOUCH eplf_gamma + call run() +end + +subroutine run + implicit none + point(1) = 0. + point(2) = 0. + integer :: i + do i=- 60,40 + point(3) = real(i)/10. + TOUCH point + print *, point(3), eplf_value, eplf_up_up, eplf_up_dn + enddo +end