From cd0cc7b1a50f15f679c3ec4091e9c145d3b4f6b2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 14 Feb 2011 12:04:15 +0100 Subject: [PATCH] Introduced second order density matrix --- src/det.irp.f | 255 ++++++++++++++++++++++++++++++++++++++++ src/eplf_function.irp.f | 224 ++--------------------------------- 2 files changed, 262 insertions(+), 217 deletions(-) diff --git a/src/det.irp.f b/src/det.irp.f index 0276d4f..6634f3e 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -255,3 +255,258 @@ BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ] END_PROVIDER +BEGIN_PROVIDER [ integer, two_e_density_num_max ] + implicit none + BEGIN_DOC +! Number of factors containing the Slater rules + END_DOC + + two_e_density_num_max = 2*mo_num + + integer :: k,l + integer :: exc(3), nact, nact2, p, p2 + integer :: det_exc + do k=1,det_num + do l=k,det_num + exc(1) = abs(det_exc(k,l,1)) + exc(2) = abs(det_exc(k,l,2)) + exc(3) = exc(1)+exc(2) + + do p=1,2 + p2 = 1+mod(p,2) + nact = elec_num_2(p) -mo_closed_num + nact2 = elec_num_2(p2)-mo_closed_num + if ( exc(3) == 0 ) then + two_e_density_num_max += 2*nact*mo_num + else if ( (exc(3) == 1).and.(exc(p) == 1) ) then + two_e_density_num_max += 2*mo_num + else if ( (exc(3) == 2).and.(exc(p) == 2) ) then + two_e_density_num_max += 2 + else if ( (exc(3) == 2).and.(exc(p) == 1) ) then + two_e_density_num_max += 1 + endif + enddo + + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, two_e_density_indice, (4,two_e_density_num_max) ] +&BEGIN_PROVIDER [ real, two_e_density_value, (2,two_e_density_num_max) ] +&BEGIN_PROVIDER [ integer, two_e_density_num ] + implicit none + BEGIN_DOC +! Compact representation of eplf factors + END_DOC + + integer :: i,j,k,l,m + + integer :: n,p,p2,q + integer :: ik,il,jk,jl, idx(4) + real :: phase + integer :: exc(4), nact, nact2 + real :: det_kl + integer :: det_exc + + two_e_density_num = 0 + + PROVIDE det + + do k=1,det_num + do l=k,det_num + + exc(1) = det_exc(k,l,1) + exc(2) = det_exc(k,l,2) + exc(4) = exc(1)*exc(2) + exc(1) = abs(exc(1)) + exc(2) = abs(exc(2)) + exc(3) = exc(1)+exc(2) + if (exc(4) /= 0) then + exc(4) = exc(4)/abs(exc(4)) + else + exc(4) = 1 + endif + phase = dble(exc(4)) + + det_kl = phase*det_coef(k)*det_coef(l) + if (k /= l) then + det_kl += det_kl + endif + + logical :: notfound +BEGIN_SHELL [ /usr/bin/python ] +code = """ + notfound = .True. + idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /) + do q=1,two_e_density_num + if (sum(abs(two_e_density_indice(:,q)-idx))) then + two_e_density_value(1,q) += det_kl + two_e_density_value(2,q) += det_kl + notfound = .False. + exit + endif + enddo + if (notfound) then + two_e_density_num += 1 + two_e_density_indice(:,two_e_density_num)=idx + two_e_density_value(1,two_e_density_num) = det_kl + two_e_density_value(2,two_e_density_num) = det_kl + endif + + notfound = .True. + idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /) + do q=1,two_e_density_num + if (sum(abs(two_e_density_indice(:,q)-idx))) then + two_e_density_value(1,q) -= det_kl + notfound = .False. + exit + endif + enddo + if (notfound) then + two_e_density_num += 1 + two_e_density_indice(:,two_e_density_num)=idx + two_e_density_value(1,two_e_density_num) = -det_kl + two_e_density_value(2,two_e_density_num) = 0. + endif +""" + +code1 = """ + idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /) + notfound = .True. + do q=1,two_e_density_num + if (sum(abs(two_e_density_indice(:,q)-idx))) then + two_e_density_value(1,q) += det_kl + notfound = .False. + exit + endif + enddo + if (notfound) then + two_e_density_num += 1 + two_e_density_indice(:,two_e_density_num)=idx + two_e_density_value(1,two_e_density_num) = det_kl + two_e_density_value(2,two_e_density_num) = 0. + endif + + notfound = .True. + idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /) + do q=1,two_e_density_num + if (sum(abs(two_e_density_indice(:,q)-idx))) then + two_e_density_value(1,q) -= det_kl + notfound = .False. + exit + endif + enddo + if (notfound) then + two_e_density_num += 1 + two_e_density_indice(:,two_e_density_num)=idx + two_e_density_value(1,two_e_density_num) = -det_kl + two_e_density_value(2,two_e_density_num) = 0. + endif +""" + +code2 = """ + notfound = .True. + idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /) + do q=1,two_e_density_num + if (sum(abs(two_e_density_indice(:,q)-idx))) then + two_e_density_value(2,q) += det_kl + notfound = .False. + exit + endif + enddo + if (notfound) then + two_e_density_num += 1 + two_e_density_indice(:,two_e_density_num)=idx + two_e_density_value(1,two_e_density_num) = 0. + two_e_density_value(2,two_e_density_num) = det_kl + endif +""" + +rep = { \ +'CLOSED' : code%{ 'I':'ik', 'J':'il', 'K':'j', 'L':'j' }, +'OPEN_CLOSED' : code%{ 'I':'j', 'J':'j', 'K':'ik', 'L':'il' }, +'OPEN_OPEN_1' : code1%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' }, +'OPEN_OPEN_2' : code2%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' } +} + + +print """ + + do p=1,2 + p2 = 1+mod(p,2) + nact = elec_num_2(p) -mo_closed_num + nact2 = elec_num_2(p2)-mo_closed_num + + if ( exc(3) == 0 ) then + do n=1,nact + ik = det(n,k,p) + il = det(n,l,p) + do j=1,mo_closed_num + ! Closed-open shell interactions + %(CLOSED)s + !- Open-closed shell interactions + %(OPEN_CLOSED)s + enddo + + !- Open-open shell interactions + do m=1,nact + jk = det(m,k,p) + jl = det(m,l,p) + %(OPEN_OPEN_1)s + enddo + do m=1,nact2 + jk = det(m,k,p2) + jl = det(m,l,p2) + %(OPEN_OPEN_2)s + 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) + + do j=1,mo_closed_num + !- Open-closed shell interactions + %(CLOSED)s + !- Closed-open shell interactions + %(OPEN_CLOSED)s + enddo + + !- Open-open shell interactions + do m=1,nact + jk = det(m,k,p) + jl = det(m,l,p) + %(OPEN_OPEN_1)s + enddo + do m=1,nact2 + jk = det(m,k,p2) + jl = det(m,l,p2) + %(OPEN_OPEN_2)s + 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) + %(OPEN_OPEN_1)s + + 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) + %(OPEN_OPEN_2)s + + endif + + enddo +"""%(rep) +END_SHELL + + enddo + enddo + +END_PROVIDER + diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 28c9bac..83c7894 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -96,228 +96,18 @@ END_PROVIDER integer :: k,l,m - do m=1,eplf_factor_num_max - i=eplf_factor_indice(1,m) - j=eplf_factor_indice(2,m) - k=eplf_factor_indice(3,m) - l=eplf_factor_indice(4,m) + do m=1,two_e_density_num + i=two_e_density_indice(1,m) + j=two_e_density_indice(2,m) + k=two_e_density_indice(3,m) + l=two_e_density_indice(4,m) temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l) - eplf_up_up += eplf_factor_value(1,m)*temp - eplf_up_dn += eplf_factor_value(2,m)*temp + eplf_up_up += two_e_density_value(1,m)*temp + eplf_up_dn += two_e_density_value(2,m)*temp enddo END_PROVIDER -BEGIN_PROVIDER [ integer, eplf_factor_num_max ] - implicit none - BEGIN_DOC -! Number of factors containing the Slater rules - END_DOC - - eplf_factor_num_max = 0 - - integer :: k,l - integer :: exc(3), nact, nact2, p, p2 - integer :: det_exc - do k=1,det_num - do l=k,det_num - exc(1) = det_exc(k,l,1) - exc(2) = det_exc(k,l,2) - exc(4) = exc(1)*exc(2) - exc(1) = abs(exc(1)) - exc(2) = abs(exc(2)) - exc(3) = exc(1)+exc(2) - - do p=1,2 - p2 = 1+mod(p,2) - nact = elec_num_2(p) -mo_closed_num - nact2 = elec_num_2(p2)-mo_closed_num - if ( exc(3) == 0 ) then - eplf_factor_num_max += 2*nact*mo_num - else if ( (exc(3) == 1).and.(exc(p) == 1) ) then - eplf_factor_num_max += 2*mo_num - else if ( (exc(3) == 2).and.(exc(p) == 2) ) then - eplf_factor_num_max += 2 - else if ( (exc(3) == 2).and.(exc(p) == 1) ) then - eplf_factor_num_max += 1 - endif - enddo - - enddo - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ integer, eplf_factor_indice, (4,eplf_factor_num_max) ] -&BEGIN_PROVIDER [ real, eplf_factor_value, (2,eplf_factor_num_max) ] - implicit none - BEGIN_DOC -! Compact representation of eplf factors - END_DOC - - integer :: i,j,k,l,m - - m=1 - do i=1,mo_num - do j=1,mo_num - do k=1,mo_num - do l=1,mo_num - if ( (eplf_factor(1,l,k,j,i) /= 0.).or. & - (eplf_factor(2,l,k,j,i) /= 0.) ) then - eplf_factor_indice(1,m) = l - eplf_factor_indice(2,m) = k - eplf_factor_indice(3,m) = j - eplf_factor_indice(4,m) = i - eplf_factor_value(1,m) = eplf_factor(1,l,k,j,i) - eplf_factor_value(2,m) = eplf_factor(2,l,k,j,i) - m += 1 - endif - enddo - enddo - enddo - enddo - FREE eplf_factor - -END_PROVIDER - -BEGIN_PROVIDER [ real, eplf_factor, (2,mo_num,mo_num,mo_num,mo_num) ] - implicit none - BEGIN_DOC -! Factors containing the Slater rules - END_DOC - - integer :: i, j - - integer :: k,l,m,n,p,p2 - integer :: ik,il,jk,jl - real :: phase - integer :: exc(4), nact, nact2 - real :: det_kl - integer :: det_exc - - do m=1,2 - do i=1,mo_num - do j=1,mo_num - do k=1,mo_num - do l=1,mo_num - eplf_factor(m,l,k,j,i) = 0. - enddo - enddo - enddo - enddo - enddo - - PROVIDE det - - do k=1,det_num - do l=k,det_num - - exc(1) = det_exc(k,l,1) - exc(2) = det_exc(k,l,2) - exc(4) = exc(1)*exc(2) - exc(1) = abs(exc(1)) - exc(2) = abs(exc(2)) - exc(3) = exc(1)+exc(2) - if (exc(4) /= 0) then - exc(4) = exc(4)/abs(exc(4)) - else - exc(4) = 1 - endif - phase = dble(exc(4)) - - det_kl = phase*det_coef(k)*det_coef(l) - if (k /= l) then - det_kl += det_kl - endif - - do p=1,2 - p2 = 1+mod(p,2) - nact = elec_num_2(p) -mo_closed_num - nact2 = elec_num_2(p2)-mo_closed_num - if ( exc(3) == 0 ) then - do n=1,nact - jk = det(n,k,p) - jl = det(n,l,p) - - do i=1,mo_closed_num - ! Closed-open shell interactions - eplf_factor(1,jk,jl,i,i) += det_kl - eplf_factor(2,jk,jl,i,i) += det_kl - eplf_factor(1,i,jl,jk,i) -= det_kl - - !- Open-closed shell interactions - eplf_factor(1,i,i,jk,jl) += det_kl - eplf_factor(2,i,i,jk,jl) += det_kl - eplf_factor(1,jk,i,i,jl) -= det_kl - enddo - - !- Open-open shell interactions - do m=1,nact - ik = det(m,k,p) - il = det(m,l,p) - eplf_factor(1,ik,il,jk,jl) += det_kl - eplf_factor(1,jk,il,ik,jl) -= det_kl - enddo - do m=1,nact2 - ik = det(m,k,p2) - il = det(m,l,p2) - eplf_factor(2,ik,il,jk,jl) += det_kl - 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) - - do i=1,mo_closed_num - !- Open-closed shell interactions - eplf_factor(1,ik,il,i,i) += det_kl - eplf_factor(2,ik,il,i,i) += det_kl - eplf_factor(1,i,il,ik,i) -= det_kl - - !- Closed-open shell interactions - eplf_factor(1,i,i,jk,jl) += det_kl - eplf_factor(2,i,i,jk,jl) += det_kl - eplf_factor(1,jk,i,i,jl) -= det_kl - enddo - - !- Open-open shell interactions - do m=1,nact - jk = det(m,k,p) - jl = det(m,l,p) - eplf_factor(1,ik,il,jk,jl) += det_kl - eplf_factor(1,jk,il,ik,jl) -= det_kl - enddo - do m=1,nact2 - jk = det(m,k,p2) - jl = det(m,l,p2) - eplf_factor(2,ik,il,jk,jl) += det_kl - 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) - eplf_factor(1,ik,il,jk,jl) += det_kl - eplf_factor(1,jk,il,ik,jl) -= det_kl - - 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) - eplf_factor(2,ik,il,jk,jl) += det_kl - - endif - enddo - - enddo - enddo - -END_PROVIDER - BEGIN_PROVIDER [ real, eplf_value_p ] implicit none