From e2ba2c91abb05275e2bcda9ba84aa3f10df03660 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 9 Jun 2010 15:10:14 +0200 Subject: [PATCH] EPLF MCSCF tested => OK ! --- bin/to_ezfio.py | 4 ++-- src/density.irp.f | 30 ++++++++++++++++++++++++------ src/det.irp.f | 18 +++++------------- src/eplf_function.irp.f | 28 ++++++++++++++++++---------- 4 files changed, 49 insertions(+), 31 deletions(-) diff --git a/bin/to_ezfio.py b/bin/to_ezfio.py index 162b16b..e2e3466 100755 --- a/bin/to_ezfio.py +++ b/bin/to_ezfio.py @@ -140,7 +140,7 @@ def write_ezfioFile(res,filename): ezfio.mo_basis_mo_coef = MoMatrix # Determinants - det_thr = 1.e-16 + det_thr = 1.e-6 closed_mos = res.closed_mos nactive = ezfio.get_mo_basis_mo_active_num() dets_a = [] @@ -166,7 +166,7 @@ def write_ezfioFile(res,filename): if len(dets_a[0]) > 0: for i in xrange(len(coef),0,-1): i -= 1 - if coef[i] == 0.: + if abs(coef[i]) < det_thr : dets_a.pop(i) dets_b.pop(i) coef.pop(i) diff --git a/src/density.irp.f b/src/density.irp.f index 57e32ba..86f5669 100644 --- a/src/density.irp.f +++ b/src/density.irp.f @@ -41,19 +41,28 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ] integer :: k,j,l, ik, il real :: buffer real :: phase + integer :: exc(4) 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 + exc(1) = abs(det_exc(k,l,1)) + exc(2) = abs(det_exc(k,l,2)) + exc(3) = exc(1)+exc(2) + exc(4) = exc(1)*exc(2) + if (exc(4) /= 0) then + exc(4) = exc(4)/abs(exc(4)) + endif + + phase = dble(exc(4)) + if (exc(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 += phase*det_coef(k)*det_coef(l)*buffer - else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,1) == 1) ) then + else if ( (exc(3) == 1).and.(exc(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 += phase*det_coef(k)*det_coef(l)*buffer @@ -80,18 +89,27 @@ BEGIN_PROVIDER [ real, density_beta_value_p ] integer :: k,j,l, ik, il real :: buffer real :: phase + integer :: exc(4) 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 + exc(1) = abs(det_exc(k,l,1)) + exc(2) = abs(det_exc(k,l,2)) + exc(3) = exc(1)+exc(2) + exc(4) = exc(1)*exc(2) + if (exc(4) /= 0) then + exc(4) = exc(4)/abs(exc(4)) + endif + + phase = dble(exc(4)) + if (exc(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 += phase*det_coef(k)*det_coef(l)*buffer - else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,2) == 1) ) then + else if ( (exc(3) == 1).and.(exc(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 += phase*det_coef(k)*det_coef(l)*buffer diff --git a/src/det.irp.f b/src/det.irp.f index f4a9979..7248cc8 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -33,10 +33,11 @@ BEGIN_PROVIDER [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ] END_PROVIDER -BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ] +BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 2) ] implicit none BEGIN_DOC -! Degree of excitation between two determinants. Indices are alpha, beta, alpha+beta, phase +! Degree of excitation between two determinants. Indices are alpha, beta +! The sign is the phase factor END_DOC integer :: p @@ -71,17 +72,9 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ] 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 l=1,det_num - det_exc(l,l,4) = 1 do k=l+1,det_num integer :: nperm nperm = 0 @@ -106,13 +99,12 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ] nperm += 1 endif enddo + det_exc(k,l,p) *= (1-2*mod( nperm, 2 )) enddo - det_exc(k,l,4) = 1-2*mod( nperm, 2 ) enddo enddo - - do p=1,4 + do p=1,2 do l=1,det_num do k=1,l-1 det_exc(k,l,p) = det_exc(l,k,p) diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 81cc9f8..6320a3f 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -4,10 +4,10 @@ BEGIN_PROVIDER [ real, eplf_gamma ] ! Value of the gaussian for the EPLF END_DOC include 'constants.F' - real :: eps + real :: eps, N + N = 0.1 eps = -real(dlog(tiny(1.d0))) - real :: N - eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps + eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps !eplf_gamma = 1.e10 !eplf_gamma = 1.e5 END_PROVIDER @@ -94,7 +94,7 @@ END_PROVIDER integer :: k,l,m,n,p,p2 integer :: ik,il,jk,jl double precision :: phase,dtemp(2) - integer :: exc + integer :: exc(4) PROVIDE det PROVIDE elec_num_2 @@ -102,13 +102,21 @@ END_PROVIDER do k=1,det_num do l=1,det_num - exc = det_exc(k,l,3) + 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 == 0 ) then + 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 @@ -151,7 +159,7 @@ END_PROVIDER enddo enddo - else if ( (exc == 1).and.(det_exc(k,l,p) == 1) ) then + 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) @@ -186,7 +194,7 @@ END_PROVIDER 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 + 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) @@ -195,7 +203,7 @@ END_PROVIDER 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 + 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) @@ -206,7 +214,7 @@ END_PROVIDER endif enddo - phase = dble(det_exc(k,l,4)) + 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)