From e0ce68a4b4b7eae51673c819cc13345acad60e53 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 4 Jun 2010 15:24:54 +0200 Subject: [PATCH] Slater rules seem to be OK. testing... --- src/density.irp.f | 12 ++- src/det.irp.f | 55 +++++++---- src/eplf_function.irp.f | 212 ++++++++++++++++++++++------------------ src/overlap.irp.f | 14 +-- src/test_1d.irp.f | 4 +- 5 files changed, 166 insertions(+), 131 deletions(-) diff --git a/src/density.irp.f b/src/density.irp.f index 5cac985..57e32ba 100644 --- a/src/density.irp.f +++ b/src/density.irp.f @@ -40,21 +40,23 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ] ! TODO vectorization integer :: k,j,l, ik, il real :: buffer + real :: phase PROVIDE det PROVIDE elec_alpha_num do k=1,det_num do l=1,det_num + phase = dble(det_exc(k,l,4)) if (det_exc(k,l,3) == 0) then buffer = 0. do i=1,elec_alpha_num-mo_closed_num buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1)) enddo - density_alpha_value_p += det_coef(k)*det_coef(l)*buffer + density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,1) == 1) ) then call get_single_excitation(k,l,ik,il,1) buffer = mo_value_p(ik)*mo_value_p(il) - density_alpha_value_p += det_coef(k)*det_coef(l)*buffer + density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer endif enddo @@ -77,20 +79,22 @@ BEGIN_PROVIDER [ real, density_beta_value_p ] ! TODO vectorization integer :: k,j,l, ik, il real :: buffer + real :: phase PROVIDE det PROVIDE elec_beta_num do k=1,det_num do l=1,det_num + phase = dble(det_exc(k,l,4)) if (det_exc(k,l,3) == 0) then buffer = 0. do i=1,elec_beta_num-mo_closed_num buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2)) enddo - density_beta_value_p += det_coef(k)*det_coef(l)*buffer + density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,2) == 1) ) then call get_single_excitation(k,l,ik,il,2) buffer = mo_value_p(ik)*mo_value_p(il) - density_beta_value_p += det_coef(k)*det_coef(l)*buffer + density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer endif enddo enddo diff --git a/src/det.irp.f b/src/det.irp.f index 422bfbb..f4a9979 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -33,10 +33,10 @@ BEGIN_PROVIDER [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ] END_PROVIDER -BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ] +BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ] implicit none BEGIN_DOC -! Degree of excitation between two determinants. The sign is the phase. +! Degree of excitation between two determinants. Indices are alpha, beta, alpha+beta, phase END_DOC integer :: p @@ -66,42 +66,57 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ] enddo - det_exc(l,k,p) = det_exc(k,l,p) enddo enddo enddo do l=1,det_num + det_exc(l,l,3) = 0 do k=l+1,det_num det_exc(k,l,3) = det_exc(k,l,1) + det_exc(k,l,2) enddo enddo ! Phase - do p=1,2 - do i=mo_closed_num,mo_num - integer :: det_pos(det_num) - do k=1,det_num - det_pos(k) = 0 - do j=1,elec_num_2(p)-mo_closed_num - if (det(j,k,p) == i) then - det_pos(k) = j + + do l=1,det_num + det_exc(l,l,4) = 1 + do k=l+1,det_num + integer :: nperm + nperm = 0 + do p=1,2 + integer :: buffer(0:mo_num-mo_closed_num) + do i=1,elec_num_2(p)-mo_closed_num + buffer(i) = det(i,k,p) + enddo + do i=1,elec_num_2(p)-mo_closed_num + if (buffer(i) /= det(i,l,p)) then + integer :: m + m=elec_num_2(p)-mo_closed_num + do j=i+1,elec_num_2(p)-mo_closed_num + if (buffer(i) == det(j,l,p)) then ! found + m=j + exit + endif + enddo + buffer(0) = buffer(i) + buffer(i) = det(m,l,p) + buffer(m) = buffer(0) + nperm += 1 endif enddo enddo - do k=1,det_num - do l=k+1,det_num - det_exc(k,l,3) *= -2*mod( (det_pos(k)+det_pos(l)), 2 )+1 - enddo - enddo + det_exc(k,l,4) = 1-2*mod( nperm, 2 ) enddo - enddo - do l=1,det_num - do k=l+1,det_num - det_exc(l,k,3) = det_exc(k,l,3) + + do p=1,4 + do l=1,det_num + do k=1,l-1 + det_exc(k,l,p) = det_exc(l,k,p) + enddo enddo enddo diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 5f828b1..81cc9f8 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -6,10 +6,8 @@ BEGIN_PROVIDER [ real, eplf_gamma ] include 'constants.F' real :: eps eps = -real(dlog(tiny(1.d0))) -!real :: N -!N = 0.1 -!eplf_gamma = (4./(3.*N)*pi*density_value_p)**(2./3.) * eps - eplf_gamma = density_value_p * eps + real :: N + eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps !eplf_gamma = 1.e10 !eplf_gamma = 1.e5 END_PROVIDER @@ -32,7 +30,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] implicit none BEGIN_DOC -! Array of all the for EPLF +! Array of all the for EPLF END_DOC integer :: i, j, k, l double precision :: t @@ -47,12 +45,15 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] do k=1,ao_num if (mo_coef(k,i) /= 0.) then do l=1,ao_num - t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k) - do j=i,mo_num - mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l) - enddo + if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then + t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k) + do j=i,mo_num + mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l) + enddo + endif enddo endif + enddo enddo @@ -78,20 +79,21 @@ END_PROVIDER PROVIDE mo_coef_transp - do j=1,mo_closed_num - do i=1,mo_closed_num - eplf_up_up += 2.d0* mo_value_p(i)* ( & + do i=1,mo_closed_num + do j=1,mo_closed_num + eplf_up_up += mo_value_p(i)* ( & mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & - mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) - eplf_up_dn += 2.d0* mo_value_p(i)*mo_value_p(i)* & + mo_value_p(j)*mo_eplf_integral_matrix(j,i) ) + eplf_up_dn += mo_value_p(i)*mo_value_p(i)* & mo_eplf_integral_matrix(j,j) enddo enddo + eplf_up_up *= 2.d0 + eplf_up_dn *= 2.d0 - integer :: k,l,m,n,p + integer :: k,l,m,n,p,p2 integer :: ik,il,jk,jl - double precision :: ckl - double precision :: phase + double precision :: phase,dtemp(2) integer :: exc PROVIDE det @@ -100,101 +102,117 @@ END_PROVIDER do k=1,det_num do l=1,det_num - ckl = det_coef(k)*det_coef(l) - exc = det_exc(k,l,3) - if ( exc < 0 ) then - phase = -1.0d0 - exc = -exc - else - phase = 1.0d0 - endif + dtemp(1) = 0.d0 + dtemp(2) = 0.d0 + do p=1,2 + p2 = 1+mod(p,2) + if ( exc == 0 ) then + ! Closed-open shell interactions + do j=1,mo_closed_num + do n=1,elec_num_2(p)-mo_closed_num + ik = det(n,k,p) + il = det(n,l,p) + dtemp(1) += mo_value_p(ik)* ( & + mo_value_p(il)*mo_eplf_integral_matrix(j,j) - & + mo_value_p(j)*mo_eplf_integral_matrix(j,il) ) + dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j) + enddo + enddo - if ( exc == 0 ) then - ! Sum over all alpha-alpha and beta-beta interactions - do p=1,2 + !- Open-closed shell interactions + do m=1,elec_num_2(p)-mo_closed_num + jk = det(m,k,p) + jl = det(m,l,p) + do i=1,mo_closed_num + dtemp(1) += mo_value_p(i)* ( & + mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - & + mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) ) + dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) + enddo + enddo + + !- Open-open shell interactions do m=1,elec_num_2(p)-mo_closed_num jk = det(m,k,p) jl = det(m,l,p) do n=1,elec_num_2(p)-mo_closed_num ik = det(n,k,p) il = det(n,l,p) - eplf_up_up += phase*ckl*mo_value_p(ik)* ( & + dtemp(1) += mo_value_p(ik)* ( & mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - & mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) ) enddo + do n=1,elec_num_2(p2)-mo_closed_num + ik = det(n,k,p2) + il = det(n,l,p2) + dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) + enddo enddo - enddo + + else if ( (exc == 1).and.(det_exc(k,l,p) == 1) ) then - ! Sum over all alpha-beta interactions - do m=1,elec_beta_num-mo_closed_num - jk = det(m,k,2) - jl = det(m,l,2) - do n=1,elec_alpha_num-mo_closed_num - ik = det(n,k,1) - il = det(n,l,1) - eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) & - + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) ) + ! Sum over only the sigma-sigma interactions involving the excitation + call get_single_excitation(k,l,ik,il,p) + + !- Open-closed shell interactions + do j=1,mo_closed_num + dtemp(1) += mo_value_p(ik)* ( & + mo_value_p(il)*mo_eplf_integral_matrix(j,j) - & + mo_value_p(j)*mo_eplf_integral_matrix(j,il) ) + dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j) enddo - enddo - else if ( exc == 1 ) then + !- Closed-open shell interactions + do i=1,mo_closed_num + dtemp(1) += mo_value_p(i)* ( & + mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - & + mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) ) + dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) + enddo + + !- Open-open shell interactions + do m=1,elec_num_2(p)-mo_closed_num + jk = det(m,k,p) + jl = det(m,l,p) + dtemp(1) += mo_value_p(ik)* ( & + mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - & + mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) ) + enddo + do m=1,elec_num_2(p2)-mo_closed_num + jk = det(m,k,p2) + jl = det(m,l,p2) + dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) + enddo + + else if ( (exc == 2).and.(det_exc(k,l,p) == 2) ) then + + ! Consider only the double excitations of same-spin electrons + call get_double_excitation(k,l,ik,il,jk,jl,p) + + dtemp(1) += mo_value_p(ik)* ( & + mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - & + mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) ) + + else if ( (exc == 2).and.(det_exc(k,l,p) == 1) ) then + + ! Consider only the double excitations of opposite-spin electrons + call get_single_excitation(k,l,ik,il,p) + call get_single_excitation(k,l,jk,jl,p2) + + dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - do p=1,2 - if ( det_exc(k,l,p) == 1 ) then - ! Sum over only the sigma-sigma interactions involving the excitation - call get_single_excitation(k,l,ik,il,p) - do m=1,elec_num_2(p)-mo_closed_num - jk = det(m,k,p) - jl = det(m,l,p) - eplf_up_up += phase*ckl*mo_value_p(ik)* ( & - mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - & - mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) ) - enddo - - ! Sum over only the sigma-(sigma_bar) interactions involving the excitation - integer :: p2 - p2 = 1+mod(p,2) - do m=1,elec_num_2(p2)-mo_closed_num - jk = det(m,k,p2) - jl = det(m,l,p2) - eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) & - + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) ) - enddo - endif - enddo - - else if (exc == 2) then - - if ( ( det_exc(k,l,1) == 2 ).or.( det_exc(k,l,2) == 2 ) ) then - - ! Consider only the double excitations of same-spin electrons - if ( det_exc(k,l,1) == 2 ) then - call get_double_excitation(k,l,ik,jk,il,jl,1) - else if ( det_exc(k,l,2) == 2 ) then - call get_double_excitation(k,l,ik,jk,il,jl,2) - endif - - eplf_up_up += phase*ckl*mo_value_p(ik)* ( & - mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - & - mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) ) - - else if ( det_exc(k,l,1) == 1 ) then - ! Consider only the double excitations of opposite-spin electrons - call get_single_excitation(k,l,ik,jk,1) - call get_single_excitation(k,l,il,jl,2) - - eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) & - + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) ) endif + enddo - endif + phase = dble(det_exc(k,l,4)) + eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1) + eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2) enddo enddo - END_PROVIDER @@ -304,7 +322,7 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center) end function - double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) +double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) implicit none include 'constants.F' @@ -348,24 +366,26 @@ end function di(ii) = 0.5d0*inv_p(2)*dble(ii) enddo - S(1,0) = (xp-xa) * S(0,0) + xab(1) = xp-xa + xab(2) = xp-xb + S(1,0) = xab(1) * S(0,0) if (i>1) then do ii=1,i-1 - S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0) + S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0) enddo endif - S(0,1) = (xp-xb) * S(0,0) + S(0,1) = xab(2) * S(0,0) if (j>1) then do jj=1,j-1 - S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1) + S(0,jj+1) = xab(2) * 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) + S(1,jj) = xab(1) * 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) + S(ii,jj) = xab(1) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1) enddo enddo diff --git a/src/overlap.irp.f b/src/overlap.irp.f index fe8ab0f..ff62a8b 100644 --- a/src/overlap.irp.f +++ b/src/overlap.irp.f @@ -8,13 +8,9 @@ BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ] integer :: i, j double precision :: ao_overlap do j=1,ao_num - do i=j,ao_num + do i=1,j ao_overlap_matrix(i,j) = ao_overlap(i,j) - enddo - enddo - do j=1,ao_num - do i=1,j-1 - ao_overlap_matrix(i,j) = ao_overlap(j,i) + ao_overlap_matrix(j,i) = ao_overlap_matrix(i,j) enddo enddo END_PROVIDER @@ -237,14 +233,14 @@ double precision function ao_overlap(i,j) 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) + integral(p,q) *= ao_coef(p,i)*ao_coef(q,j) enddo enddo - ao_overlap = 0. + ao_overlap = 0.d0 do q=1,ao_prim_num(j) do p=1,ao_prim_num(i) - ao_overlap = ao_overlap + integral(p,q) + ao_overlap += integral(p,q) enddo enddo diff --git a/src/test_1d.irp.f b/src/test_1d.irp.f index 37d05a1..3873c92 100644 --- a/src/test_1d.irp.f +++ b/src/test_1d.irp.f @@ -10,8 +10,8 @@ subroutine run point(1) = 0. point(2) = 0. integer :: i - do i=0,40 - point(3) = real(i)/10. + do i=-40,40 + point(3) = real(i)/20. TOUCH point print *, point(3), eplf_value_p, eplf_up_up, eplf_up_dn enddo