fixed a bug in two_rdm, added the possibility to Write/Read the all_states active 2 rdm

src/bitmask/EZFIO.cfg

@@ -0,0 +1,5 @@
+[n_act_orb] 
+type: integer
+doc: Number of active |MOs| 
+interface: ezfio
+

src/bitmask/core_inact_act_virt.irp.f

@@ -49,9 +49,10 @@ BEGIN_PROVIDER [ integer, n_act_orb]
       n_act_orb += 1
     endif
   enddo
-  
   call write_int(6,n_act_orb, 'Number of active   MOs')
-   
+  if (mpi_master) then
+    call ezfio_set_bitmask_n_act_orb(n_act_orb)
+  endif
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer, n_virt_orb ]

src/two_body_rdm/EZFIO.cfg

@@ -0,0 +1,48 @@
+[two_rdm_ab_disk]  
+type: double precision
+doc: active part of the two body rdm alpha/beta stored on disk
+interface: ezfio
+size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
+
+[io_two_body_rdm_ab]
+type: Disk_access
+doc: Read/Write the active part of the two-body rdm for alpha/beta electrons from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
+[two_rdm_aa_disk]  
+type: double precision
+doc: active part of the two body rdm alpha/alpha stored on disk
+interface: ezfio
+size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
+
+[io_two_body_rdm_aa]
+type: Disk_access
+doc: Read/Write the active part of the two-body rdm for alpha/alpha electrons from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
+[two_rdm_bb_disk]  
+type: double precision
+doc: active part of the two body rdm beta/beta stored on disk
+interface: ezfio
+size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
+
+[io_two_body_rdm_bb]
+type: Disk_access
+doc: Read/Write the active part of the two-body rdm for beta/beta electrons from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
+[two_rdm_spin_trace_disk]  
+type: double precision
+doc: active part of the two body rdm spin trace stored on disk
+interface: ezfio
+size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
+
+[io_two_body_rdm_spin_trace]
+type: Disk_access
+doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+

src/two_body_rdm/README.rst

@@ -3,7 +3,6 @@ two_body_rdm
 ============
 
 Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
-arrays, with pysicists notation, consistent with the two-electron integrals in the
+arrays, with pysicists notation, consistent with the two-electron integrals in the MO basis.
-MO basis.
 
 

src/two_body_rdm/act_2_rdm.irp.f

@@ -26,8 +26,17 @@
  ispin = 3 
  act_2_rdm_ab_mo = 0.d0
  call wall_time(wall_1)
- call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ if(read_two_body_rdm_ab)then
-
+  print*,'Reading act_2_rdm_ab_mo from disk ...'
+  call ezfio_get_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo)
+ else 
+  call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ endif
+ if(write_two_body_rdm_ab)then
+  print*,'Writing act_2_rdm_ab_mo on disk ...'
+  call ezfio_set_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo)
+  call ezfio_set_two_body_rdm_io_two_body_rdm_ab("Read")
+ endif
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1
  END_PROVIDER 
@@ -54,7 +63,17 @@
  ispin = 1 
  act_2_rdm_aa_mo = 0.d0
  call wall_time(wall_1)
- call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ if(read_two_body_rdm_aa)then
+  print*,'Reading act_2_rdm_aa_mo from disk ...'
+  call ezfio_get_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo)
+ else 
+  call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ endif
+ if(write_two_body_rdm_aa)then
+  print*,'Writing act_2_rdm_aa_mo on disk ...'
+  call ezfio_set_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo)
+  call ezfio_set_two_body_rdm_io_two_body_rdm_aa("Read")
+ endif
 
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1
@@ -82,7 +101,17 @@
  ispin = 2 
  act_2_rdm_bb_mo = 0.d0
  call wall_time(wall_1)
- call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ if(read_two_body_rdm_bb)then
+  print*,'Reading act_2_rdm_bb_mo from disk ...'
+  call ezfio_get_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo)
+ else 
+  call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ endif
+ if(write_two_body_rdm_bb)then
+  print*,'Writing act_2_rdm_bb_mo on disk ...'
+  call ezfio_set_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo)
+  call ezfio_set_two_body_rdm_io_two_body_rdm_bb("Read")
+ endif
 
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1
@@ -109,7 +138,17 @@
  ispin = 4 
  act_2_rdm_spin_trace_mo = 0.d0
  call wall_time(wall_1)
- call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ if(read_two_body_rdm_spin_trace)then
+  print*,'Reading act_2_rdm_spin_trace_mo from disk ...'
+  call ezfio_get_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo)
+ else 
+  call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ endif
+ if(write_two_body_rdm_spin_trace)then
+  print*,'Writing act_2_rdm_spin_trace_mo on disk ...'
+  call ezfio_set_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo)
+  call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
+ endif
 
  call wall_time(wall_2)
  print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1

src/two_rdm_routines/davidson_like_2rdm.irp.f

@@ -137,7 +137,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
-   sze_buff = 6 * norb
+   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 
    list_orb_reverse = -1000 
    do i = 1, norb
     list_orb_reverse(list_orb(i)) = i 
@@ -267,7 +267,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
           endif
          endif
          call orb_range_off_diag_double_to_all_states_ab_dm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-!         print*,'coucou'
+
        enddo
       endif
        
@@ -347,15 +347,13 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
          nkeys = 0
        endif
-        call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
+       call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        ! increment the alpha/alpha part for single excitations
        if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
          nkeys = 0
        endif
-!       call orb_range_off_diag_single_to_2_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
+       call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-        call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-!       print*,'to do orb_range_off_diag_single_to_2_rdm_aa_dm_buffer'
        endif
        
      enddo
@@ -382,7 +380,6 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
           nkeys = 0
         endif
         call orb_range_off_diag_double_to_all_states_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-!        print*,'to do orb_range_off_diag_double_to_2_rdm_aa_dm_buffer'
       enddo
      endif
      call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
@@ -453,7 +450,6 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
           nkeys = 0
         endif
          call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-!        print*,'to do orb_range_off_diag_single_to_2_rdm_ab_dm_buffer'
         ! increment the beta /beta  part for single excitations
         if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
           call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)

src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f

@@ -139,7 +139,7 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
    PROVIDE N_int
 
    call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
-   sze_buff = norb ** 3 + 6 * norb
+   sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 
    list_orb_reverse = -1000 
    do i = 1, norb
     list_orb_reverse(list_orb(i)) = i 
@@ -271,11 +271,12 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
           endif
          endif
          call orb_range_off_diag_double_to_2_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-!         print*,'todo orb_range_off_diag_double_to_2_rdm_ab_dm_buffer'
 
        enddo
       endif
        
+      call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+      nkeys = 0
      enddo
      
    enddo
@@ -364,6 +365,8 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
        
      enddo
      
+     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     nkeys = 0
      
      ! Compute Hij for all alpha doubles
      ! ----------------------------------
@@ -389,6 +392,8 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
         call orb_range_off_diag_double_to_2_rdm_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
       enddo
      endif
+     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     nkeys = 0
      
      
      ! Single and double beta excitations
@@ -466,6 +471,8 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
         call orb_range_off_diag_single_to_2_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        endif
      enddo
+     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     nkeys = 0
      
      ! Compute Hij for all beta doubles
      ! ----------------------------------

src/two_rdm_routines/update_rdm.irp.f

@@ -257,11 +257,9 @@
  if(list_orb_reverse(p2).lt.0)return
  p2 = list_orb_reverse(p2)
  if(alpha_beta)then
-!   print*,'coucou'
    nkeys += 1
    do istate = 1, N_st
     values(istate,nkeys) = c_1(istate) * phase
-!    print*,'values',values(istate,nkeys),nkeys
    enddo
    keys(1,nkeys) = h1
    keys(2,nkeys) = h2