Checkpoint of density matrices

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)

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
 
 

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

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 = [\