diff --git a/eplf.config b/eplf.config index 33c656f..3ee4734 100644 --- a/eplf.config +++ b/eplf.config @@ -41,6 +41,12 @@ grid num_y integer = at(grid_point_num,2) num_z integer = at(grid_point_num,3) +density_matrix + one real (mo_basis_mo_active_num,mo_basis_mo_active_num,2) + two_num integer + two_indice real (4,density_matrix_two_num) + two_value real (2,density_matrix_two_num) + grid_data eplf real (grid_num_x,grid_num_y,grid_num_z) eplf_grad real (grid_num_x,grid_num_y,grid_num_z,4) diff --git a/src/density.irp.f b/src/density.irp.f index 394d738..4fc1d34 100644 --- a/src/density.irp.f +++ b/src/density.irp.f @@ -37,10 +37,10 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ] density_alpha_value_p += mo_value_p(i)**2 enddo - do m=1,one_e_density_num - i = one_e_density_indice(1,m) - j = one_e_density_indice(2,m) - density_alpha_value_p += mo_value_prod_p(i,j) * one_e_density_value(1,m) + do j=1,elec_alpha_num-mo_closed_num + do i=1,elec_alpha_num-mo_closed_num + density_alpha_value_p += mo_value_prod_p(i,j) * one_e_density_mo(i,j,1) + enddo enddo END_PROVIDER @@ -52,17 +52,16 @@ BEGIN_PROVIDER [ real, density_beta_value_p ] END_DOC density_beta_value_p = 0. - integer :: i + do i=1,mo_closed_num density_beta_value_p += mo_value_p(i)**2 enddo - do m=1,one_e_density_num - i = one_e_density_indice(1,m) - j = one_e_density_indice(2,m) - density_beta_value_p += mo_value_prod_p(i,j) * one_e_density_value(2,m) + do j=1,elec_beta_num-mo_closed_num + do i=1,elec_beta_num-mo_closed_num + density_beta_value_p += mo_value_prod_p(i,j) * one_e_density_mo(i,j,2) + enddo enddo - END_PROVIDER diff --git a/src/det.irp.f b/src/det.irp.f index 07d45cb..83b2662 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -55,44 +55,41 @@ integer function det_exc(k,l,p) jmax = elec_num_2(p)-mo_closed_num det_exc = 0 + integer :: dl(mo_closed_num), dk(mo_closed_num), buffer(0:mo_closed_num) do i=1,jmax - logical :: found - found = .False. + dk(i) = det(i,k,p) + dl(i) = det(i,l,p) + buffer(i) = dk(i) + enddo + + integer :: kmax + logical :: notfound + do i=1,jmax + notfound = .True. do j=1,jmax - if (det(j,l,p) == det(i,k,p)) then - found = .True. - exit - endif + notfound = notfound .and. (dl(j) /= dk(i)) enddo - if (.not.found) then + if (notfound) then det_exc += 1 endif - enddo ! Phase integer :: nperm nperm = 0 - integer, allocatable, save :: buffer(:) - if (.not. allocated(buffer)) then - allocate (buffer(0:mo_num-mo_closed_num)) - endif - 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 + do i=1,jmax + if (buffer(i) /= dl(i)) 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=jmax + do j=i+1,jmax + if (buffer(i) == dl(j)) then ! found m=j exit endif enddo buffer(0) = buffer(i) - buffer(i) = det(m,l,p) + buffer(i) = dl(m) buffer(m) = buffer(0) nperm += m-i endif @@ -108,33 +105,37 @@ subroutine get_single_excitation(k,l,m,n,p) integer, intent(out) :: m, n ! m->n excitation integer, intent(in) :: p ! spin - logical :: found + logical :: notfound integer :: i,j m=0 n=0 - do j=1,elec_num_2(p)-mo_closed_num - found = .False. - do i=1,elec_num_2(p)-mo_closed_num - if (det(j,k,p) == det(i,l,p)) then - found = .True. - exit - endif + + integer :: dl(mo_closed_num), dk(mo_closed_num), buffer(0:mo_closed_num) + integer :: jmax + jmax = elec_num_2(p)-mo_closed_num + + do i=1,jmax + dk(i) = det(i,k,p) + dl(i) = det(i,l,p) + enddo + + do j=1,jmax + notfound = .True. + do i=1,jmax + notfound = notfound .and. (dk(j) /= dl(i)) enddo - if (.not.found) then - m = det(j,k,p) + if (notfound) then + m = dk(j) exit endif enddo - do i=1,elec_num_2(p)-mo_closed_num - found = .False. - do j=1,elec_num_2(p)-mo_closed_num - if (det(j,k,p) == det(i,l,p)) then - found = .True. - exit - endif + do i=1,jmax + notfound = .True. + do j=1,jmax + notfound = notfound .and. (dk(j) /= dl(i)) enddo - if (.not.found) then + if (notfound) then n = det(i,l,p) exit endif @@ -149,43 +150,47 @@ subroutine get_double_excitation(k,l,m,n,r,s,p) integer, intent(out) :: r, s ! r->s excitation integer, intent(in) :: p ! spin - logical :: found + logical :: notfound integer :: i,j m=0 n=0 r=0 s=0 - do j=1,elec_num_2(p)-mo_closed_num - found = .False. - do i=1,elec_num_2(p)-mo_closed_num - if (det(j,k,p) == det(i,l,p)) then - found = .True. - exit - endif + + integer :: dl(mo_closed_num), dk(mo_closed_num), buffer(0:mo_closed_num) + integer :: jmax + jmax = elec_num_2(p)-mo_closed_num + + do i=1,jmax + dk(i) = det(i,k,p) + dl(i) = det(i,l,p) + enddo + + do j=1,jmax + notfound = .True. + do i=1,jmax + notfound = notfound .and. (dk(j) /= dl(i)) enddo - if (.not.found) then + if (notfound) then if (m == 0) then - m = det(j,k,p) + m = dk(j) else - r = det(j,k,p) + r = dk(j) exit endif endif enddo - do i=1,elec_num_2(p)-mo_closed_num - found = .False. - do j=1,elec_num_2(p)-mo_closed_num - if (det(j,k,p) == det(i,l,p)) then - found = .True. - exit - endif + do j=1,jmax + notfound = .True. + do i=1,jmax + notfound = notfound .and. (dk(j) /= dl(i)) enddo - if (.not.found) then + if (notfound) then if (n == 0) then - n = det(i,l,p) + n = dl(j) else - s = det(i,l,p) + s = dl(j) exit endif endif @@ -193,72 +198,6 @@ subroutine get_double_excitation(k,l,m,n,r,s,p) end -BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ] - implicit none - BEGIN_DOC -! Spin Density matrix in the AO basis - END_DOC - integer :: i,j,k,l,m,ispin, ik,il - do ispin=1,3 - - do j=1,mo_num - do i=1,mo_num - ci_mo(i,j,ispin) = 0. - enddo - enddo - - do l=1,det_num - do m=1,det_num - real :: factor - factor = 2.*det_coef(l)*det_coef(m) - do il=1,mo_closed_num - do ik=1,mo_closed_num - ci_mo(ik,il,ispin) += factor - enddo - enddo - enddo - enddo - - enddo - - - do l=1,det_num - do m=1,det_num - factor = det_coef(l)*det_coef(m) - do ispin=1,2 - do j=mo_closed_num+1,elec_num_2(ispin) - ik = det(j-mo_closed_num,l,ispin) - do il=1,mo_closed_num - ci_mo(ik,il,ispin) += factor - ci_mo(il,ik,ispin) += factor - enddo - do i=mo_closed_num+1,elec_num_2(ispin) - il = det(i-mo_closed_num,m,ispin) - ci_mo(ik,il,ispin) += factor - ci_mo(il,ik,ispin) += factor - enddo - enddo - enddo - - ispin=3 - do j=mo_closed_num+1,elec_num_2(1) - ik = det(j-mo_closed_num,l,1) - do il=1,mo_closed_num - ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) - ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) - enddo - do i=mo_closed_num+1,elec_num_2(2) - il = det(i-mo_closed_num,m,2) - ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) - ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) - enddo - enddo - - enddo - enddo - -END_PROVIDER - BEGIN_PROVIDER [ integer, two_e_density_num_max ] implicit none @@ -266,8 +205,13 @@ BEGIN_PROVIDER [ integer, two_e_density_num_max ] ! Number of factors containing the Slater rules END_DOC - two_e_density_num_max = 2*mo_num + two_e_density_num_max = 0 + call get_density_matrix_two_num(two_e_density_num_max) + if (two_e_density_num_max /= 0) then + return + endif + two_e_density_num_max = 2*mo_num integer :: k,l integer :: exc(3), nact, nact2, p, p2 integer :: det_exc @@ -295,6 +239,7 @@ BEGIN_PROVIDER [ integer, two_e_density_num_max ] enddo enddo + call set_density_matrix_two_num(two_e_density_num_max) END_PROVIDER BEGIN_PROVIDER [ integer, two_e_density_indice, (4,two_e_density_num_max) ] @@ -305,6 +250,13 @@ END_PROVIDER ! Compact representation of eplf factors END_DOC + two_e_density_indice(1,1) = -1 + call get_density_matrix_two_indice(two_e_density_indice) + call get_density_matrix_two_value(two_e_density_value) + if (two_e_density_indice(1,1) /= -1) then + return + endif + integer :: i,j,k,l,m integer :: n,p,p2,q @@ -509,6 +461,10 @@ END_SHELL enddo enddo + call set_density_matrix_two_indice(two_e_density_indice) + call set_density_matrix_two_value(two_e_density_value) + call set_density_matrix_two_num(two_e_density_num) + END_PROVIDER BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ] @@ -517,6 +473,13 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ] ! One electron spin density matrix in MO space END_DOC integer :: i,j,k,l,p, il, jl + + one_e_density_mo(1,1,1) = -1. + call get_density_matrix_one(one_e_density_mo) + if (one_e_density_mo(1,1,1) /= -1.) then + return + endif + do p=1,2 do i=1,mo_active_num do j=1,mo_active_num @@ -554,4 +517,6 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ] enddo enddo + call set_density_matrix_one(one_e_density_mo) + END_PROVIDER diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index c6305f8..05905d6 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -45,6 +45,10 @@ data = [ \ ("grid_data_elf_partition" , "real" , "(grid_x_num,grid_y_num,grid_z_num)" ), ("grid_data_eplf_partition" , "real" , "(grid_x_num,grid_y_num,grid_z_num)" ), ("grid_data_density_partition" , "real" , "(grid_x_num,grid_y_num,grid_z_num)" ), +("density_matrix_one" , "real" , "(mo_active_num,mo_active_num,2)" ), +("density_matrix_two_num" , "integer" , "" ), +("density_matrix_two_indice" , "real" , "(4,two_e_density_num_max)" ), +("density_matrix_two_value" , "real" , "(2,two_e_density_num_max)" ), ] data_no_set = [\