BEGIN_PROVIDER [ real, eplf_gamma ] implicit none BEGIN_DOC ! Value of the gaussian for the EPLF END_DOC include 'constants.F' real :: eps, N ! N = 0.1 ! eps = -real(dlog(tiny(1.d0))) ! eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps N=0.001 eps = 50. eplf_gamma = eps**2 *0.5d0*(4./9.*pi*density_value_p * N**(-1./3.))**(2./3.) !eplf_gamma = 1.e10 !eplf_gamma = 1.e5 END_PROVIDER BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ] implicit none BEGIN_DOC ! Array of all the for EPLF END_DOC integer :: i, j double precision :: ao_eplf_integral do i=1,ao_num do j=i,ao_num ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,dble(eplf_gamma),point) ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i) enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] implicit none BEGIN_DOC ! Array of all the for EPLF END_DOC integer :: i, j, k, l double precision :: t PROVIDE ao_eplf_integral_matrix PROVIDE mo_coef do i=1,mo_num do j=i,mo_num mo_eplf_integral_matrix(j,i) = 0.d0 enddo do k=1,ao_num if (mo_coef(k,i) /= 0.) then do l=1,ao_num 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 do i=1,mo_num do j=i+1,mo_num mo_eplf_integral_matrix(i,j) = mo_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 eplf_up_up = 0.d0 eplf_up_dn = 0.d0 PROVIDE mo_coef_transp 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(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,p2 integer :: ik,il,jk,jl double precision :: phase,dtemp(2) integer :: exc(4) PROVIDE det PROVIDE elec_num_2 do k=1,det_num do l=1,det_num exc(1) = abs(det_exc(k,l,1)) exc(2) = abs(det_exc(k,l,2)) exc(3) = exc(1)+exc(2) exc(4) = det_exc(k,l,1)*det_exc(k,l,2) if (exc(4) /= 0) then exc(4) = exc(4)/abs(exc(4)) else exc(4) = 1 endif dtemp(1) = 0.d0 dtemp(2) = 0.d0 do p=1,2 p2 = 1+mod(p,2) if ( exc(3) == 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 !- 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) 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 else if ( (exc(3) == 1).and.(exc(p) == 1) ) then ! 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 !- 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(3) == 2).and.(exc(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(3) == 2).and.(exc(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) endif enddo phase = dble(exc(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 BEGIN_PROVIDER [ real, eplf_value_p ] implicit none BEGIN_DOC ! Value of the EPLF at the current point. END_DOC double precision :: aa, ab double precision, parameter :: eps = tiny(1.d0) aa = eplf_up_up ab = eplf_up_dn if ( (aa > 0.d0).and.(ab > 0.d0) ) then aa = min(1.d0,aa) ab = min(1.d0,ab) aa = -dlog(aa) aa = dsqrt(aa) ab = -dlog(ab) ab = dsqrt(ab) eplf_value_p = (aa-ab)/(aa+ab+eps) else eplf_value_p = 0.d0 endif END_PROVIDER !double precision function ao_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=10000 ! 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 += & ! (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) ! x = x+dx ! enddo ! ao_eplf_integral_primitive_oneD_numeric = dtemp*dx ! !end function ! !double precision function ao_eplf_integral_numeric(i,j,gmma,center) ! implicit none ! integer, intent(in) :: i, j ! integer :: p,q,k ! double precision :: integral ! double precision :: ao_eplf_integral_primitive_oneD_numeric ! real :: gmma, center(3), c ! ! ! 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) ! integral = & ! ao_eplf_integral_primitive_oneD_numeric( & ! ao_expo(i,p), & ! nucl_coord(ao_nucl(i),1), & ! ao_power(i,1), & ! ao_expo(j,q), & ! nucl_coord(ao_nucl(j),1), & ! ao_power(j,1), & ! gmma, & ! center(1)) * & ! ao_eplf_integral_primitive_oneD_numeric( & ! ao_expo(i,p), & ! nucl_coord(ao_nucl(i),2), & ! ao_power(i,2), & ! ao_expo(j,q), & ! nucl_coord(ao_nucl(j),2), & ! ao_power(j,2), & ! gmma, & ! center(2)) * & ! ao_eplf_integral_primitive_oneD_numeric( & ! ao_expo(i,p), & ! nucl_coord(ao_nucl(i),3), & ! ao_power(i,3), & ! ao_expo(j,q), & ! nucl_coord(ao_nucl(j),3), & ! ao_power(j,3), & ! gmma, & ! center(3)) ! ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral ! enddo ! enddo ! !end function double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) implicit none include 'constants.F' real, intent(in) :: a,b ! Exponents double precision , intent(in) :: gmma ! eplf_gamma real, intent(in) :: xa,xb,xr ! Centers integer, intent(in) :: i,j ! Powers of xa and xb integer :: ii, jj, kk ASSERT (a>0.) ASSERT (b>0.) ASSERT (i>=0) ASSERT (j>=0) double precision :: alpha, beta, delta, c alpha = a*xa*xa + b*xb*xb + gmma*xr*xr delta = a*xa + b*xb + gmma*xr beta = a + b + gmma c = alpha-delta**2/beta if ( c > 32.d0 ) then ! Cut-off on exp(-32) ao_eplf_integral_primitive_oneD = 0.d0 return endif double precision :: u,v,w c = exp(-c) w = 1.d0/sqrt(beta) u=delta/beta-xa v=delta/beta-xb double precision :: fact2, binom, f, ac ac = 0.d0 integer :: istart(0:1) istart = (/ 0, 1 /) do ii=0,i do jj=istart(mod(ii,2)),j,2 kk=(ii+jj)/2 f=binom(ii,i)*binom(jj,j)*fact2(ii+jj-1) ac += u**(i-ii)*v**(j-jj)*w**(ii+jj)*f/dble(2**kk) enddo enddo ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac end function double precision function ao_eplf_integral(i,j,gmma,center) implicit none integer, intent(in) :: i, j real, intent(in) :: center(3) double precision, intent(in) :: gmma !DEC$ ATTRIBUTES FORCEINLINE integer :: p,q,k double precision :: integral double precision :: ao_eplf_integral_primitive_oneD double precision :: buffer(100) ASSERT(i>0) ASSERT(j>0) ASSERT(ao_prim_num_max < 100) ASSERT(i<=ao_num) ASSERT(j<=ao_num) ao_eplf_integral = 0.d0 do p=1,ao_prim_num_max buffer(p) = 0.d0 enddo do q=1,ao_prim_num(j) do p=1,ao_prim_num(i) integral = & ao_eplf_integral_primitive_oneD( & ao_expo(i,p), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & ao_expo(j,q), & nucl_coord(ao_nucl(j),1), & ao_power(j,1), & gmma, & center(1)) if (integral /= 0.d0) then integral *= & ao_eplf_integral_primitive_oneD( & ao_expo(i,p), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & ao_expo(j,q), & nucl_coord(ao_nucl(j),2), & ao_power(j,2), & gmma, & center(2)) if (integral /= 0.d0) then integral *= & ao_eplf_integral_primitive_oneD( & ao_expo(i,p), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & ao_expo(j,q), & nucl_coord(ao_nucl(j),3), & ao_power(j,3), & gmma, & center(3)) buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q) endif endif enddo enddo do p=1,ao_prim_num_max ao_eplf_integral += buffer(p) enddo end function double precision function mo_eplf_integral(i,j) implicit none integer :: i, j, k, l PROVIDE ao_eplf_integral_matrix PROVIDE mo_coef mo_eplf_integral = 0.d0 do k=1,ao_num if (mo_coef(k,i) /= 0.) then do l=1,ao_num mo_eplf_integral += & mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l) enddo endif enddo end function