Integrals AO not recalcaulated if not necessary

src/BiInts/map_integrals.irp.f

@@ -382,6 +382,7 @@ subroutine dump_$ao_integrals(filename)
   integer(cache_key_kind), pointer :: key(:)
   real(integral_kind), pointer   :: val(:)
   integer*8                      :: i,j, n
+  call ezfio_set_work_empty(.False.)
   open(unit=66,file=filename,FORM='unformatted')
   write(66) integral_kind, key_kind
   write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size,    &

src/BiInts/mo_bi_integrals.irp.f

@@ -319,14 +319,14 @@ end
 
 
 
- BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ]
+ BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_from_ao, (mo_tot_num_align,mo_tot_num) ]
-&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ]
-&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti_from_ao, (mo_tot_num_align,mo_tot_num) ]
   implicit none
   BEGIN_DOC
-  ! mo_bielec_integral_jj(i,j) = J_ij
+  ! mo_bielec_integral_jj_from_ao(i,j) = J_ij
-  ! mo_bielec_integral_jj_exchange(i,j) = J_ij
+  ! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij
-  ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij
+  ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij
   END_DOC
   
   integer                        :: i,j,p,q,r,s
@@ -342,8 +342,8 @@ end
     PROVIDE ao_bielec_integrals_in_map
   endif
   
-  mo_bielec_integral_jj = 0.d0
+  mo_bielec_integral_jj_from_ao = 0.d0
-  mo_bielec_integral_jj_exchange = 0.d0
+  mo_bielec_integral_jj_exchange_from_ao = 0.d0
   
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
   
@@ -353,7 +353,7 @@ end
       !$OMP  iqrs, iqsr,iqri,iqis)                                   &
       !$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num,&
       !$OMP  ao_integrals_threshold,do_direct_integrals,abort_here)  &
-      !$OMP REDUCTION(+:mo_bielec_integral_jj,mo_bielec_integral_jj_exchange)
+      !$OMP REDUCTION(+:mo_bielec_integral_jj_from_ao,mo_bielec_integral_jj_exchange_from_ao)
   
   allocate( int_value(ao_num), int_idx(ao_num),                      &
       iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),&
@@ -439,8 +439,8 @@ end
         !DIR$ VECTOR ALIGNED
         do j=1,mo_tot_num
           c = mo_coef_transp(j,q)*mo_coef_transp(j,s)
-          mo_bielec_integral_jj(j,i) += c * iqis(i)
+          mo_bielec_integral_jj_from_ao(j,i) += c * iqis(i)
-          mo_bielec_integral_jj_exchange(j,i) += c * iqri(i)
+          mo_bielec_integral_jj_exchange_from_ao(j,i) += c * iqri(i)
         enddo
       enddo
       
@@ -453,8 +453,35 @@ end
     stop 'Aborting in MO integrals calculation'
   endif
   
-  mo_bielec_integral_jj_anti = mo_bielec_integral_jj - mo_bielec_integral_jj_exchange
+  mo_bielec_integral_jj_anti_from_ao = mo_bielec_integral_jj_from_ao - mo_bielec_integral_jj_exchange_from_ao
   
   
 END_PROVIDER
 
+
+ BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ]
+  implicit none
+  BEGIN_DOC
+  ! mo_bielec_integral_jj(i,j) = J_ij
+  ! mo_bielec_integral_jj_exchange(i,j) = K_ij
+  ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: get_mo_bielec_integral
+  
+  PROVIDE mo_bielec_integrals_in_map
+
+  mo_bielec_integral_jj = 0.d0
+  mo_bielec_integral_jj_exchange = 0.d0
+  do j=1,mo_tot_num
+   do i=1,mo_tot_num
+     mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
+     mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
+     mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) 
+   enddo
+ enddo
+  
+END_PROVIDER

src/Dets/NEEDED_MODULES

@@ -1 +1 @@
-AOs BiInts Bitmask Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils
+AOs BiInts Bitmask Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

src/Dets/determinants.irp.f

@@ -352,7 +352,7 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,psicoef)
   call ezfio_set_determinants_psi_det(psi_det_save)
   deallocate (psi_det_save)
 
-  progress_bar(7) = 7
+  progress_bar(1) = 7
   progress_value = dble(progress_bar(1))
   allocate (psi_coef_save(ndet,nstates))
   do k=1,nstates

src/Dets/ref_bitmask.irp.f

@@ -1,4 +1,5 @@
- BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ]
+ BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]
+&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ]
 &BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ]
 &BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ]
 &BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ]
@@ -15,17 +16,20 @@
   call bitstring_to_list(ref_bitmask(1,2), occ(1,2), i, N_int)
   
   
+  ref_bitmask_energy = 0.d0
   mono_elec_ref_bitmask_energy = 0.d0
   kinetic_ref_bitmask_energy   = 0.d0
   nucl_elec_ref_bitmask_energy = 0.d0
   bi_elec_ref_bitmask_energy = 0.d0
   
   do i = 1, elec_beta_num
+    ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) + mo_mono_elec_integral(occ(i,2),occ(i,2))
     kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) + mo_kinetic_integral(occ(i,2),occ(i,2))
     nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) + mo_nucl_elec_integral(occ(i,2),occ(i,2))
   enddo
 
   do i = elec_beta_num+1,elec_alpha_num
+    ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1))
     kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) 
     nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) 
   enddo
@@ -33,39 +37,21 @@
   do j= 1, elec_alpha_num
     do i = j+1, elec_alpha_num
       bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1))
+      ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1))
     enddo
   enddo
   
   do j= 1, elec_beta_num
     do i = j+1, elec_beta_num
       bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2))
+      ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2))
     enddo
     do i= 1, elec_alpha_num
       bi_elec_ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2))
+      ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2))
     enddo
   enddo
   mono_elec_ref_bitmask_energy = kinetic_ref_bitmask_energy +   nucl_elec_ref_bitmask_energy
   
 END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Energy of the reference bitmask used in Slater rules
-  END_DOC
-  
-  integer                        :: i,j
-  
-  ref_bitmask_energy = 0.d0
-  
-  do j=1,ao_num
-   do i=1,ao_num
-     ref_bitmask_energy += 0.5d0 * ( &
-         (ao_mono_elec_integral(i,j) + Fock_matrix_alpha_ao(i,j) ) *  HF_density_matrix_ao_alpha(i,j) + &
-         (ao_mono_elec_integral(i,j) + Fock_matrix_beta_ao (i,j) ) *  HF_density_matrix_ao_beta (i,j) )
-   enddo
-  enddo
-  
-END_PROVIDER
-

src/Dets/save_fock_orb.irp.f

@@ -1,45 +0,0 @@
-program save_fock_orb
-
-  double precision, allocatable :: one_body_dm_ao_alpha(:,:)
-  double precision, allocatable :: one_body_dm_ao_beta (:,:)
-
-  read_wf = .True.
-  touch read_wf
-  print *,  'N_det = ',N_det
-  print *, HF_energy
-  allocate (one_body_dm_ao_alpha(ao_num_align,ao_num),               &
-            one_body_dm_ao_beta (ao_num_align,ao_num))
-
-  call mo_to_ao_no_overlap(one_body_dm_mo_alpha,mo_tot_num_align, &
-                one_body_dm_ao_alpha,ao_num_align)
-  call mo_to_ao_no_overlap(one_body_dm_mo_beta ,mo_tot_num_align, &
-                one_body_dm_ao_beta ,ao_num_align)
-
-  do i=1,mo_tot_num
-   do j=1,mo_tot_num
-     if (abs(fock_matrix_mo(i,j)) > 1.d-10) then
-       print *, i,j, fock_matrix_mo(i,j)
-     endif
-   enddo
-  enddo
-  hf_density_matrix_ao_alpha = one_body_dm_ao_alpha
-  hf_density_matrix_ao_beta  = one_body_dm_ao_beta
-  touch hf_density_matrix_ao_alpha hf_density_matrix_ao_beta
-
-  print *,  '---'
-  do i=1,mo_tot_num
-   do j=1,mo_tot_num
-     if (abs(fock_matrix_mo(i,j)) > 1.d-10) then
-       print *, i,j, fock_matrix_mo(i,j)
-     endif
-   enddo
-  enddo
-  mo_coef = eigenvectors_fock_matrix_mo
-  mo_label = 'CASSCF'
-  TOUCH mo_coef mo_label 
-  print *, HF_energy
-  call save_mos
-
-  deallocate(one_body_dm_ao_alpha,one_body_dm_ao_beta)
-end
-

src/Hartree_Fock/Fock_matrix.irp.f

@@ -272,7 +272,17 @@ BEGIN_PROVIDER [ double precision, HF_energy ]
  BEGIN_DOC
  ! Hartree-Fock energy
  END_DOC
- HF_energy = nuclear_repulsion + ref_bitmask_energy
+ HF_energy = nuclear_repulsion
+ 
+ integer                        :: i,j
+ do j=1,ao_num
+   do i=1,ao_num
+     HF_energy += 0.5d0 * (                                          &
+         (ao_mono_elec_integral(i,j) + Fock_matrix_alpha_ao(i,j) ) *  HF_density_matrix_ao_alpha(i,j) +&
+         (ao_mono_elec_integral(i,j) + Fock_matrix_beta_ao (i,j) ) *  HF_density_matrix_ao_beta (i,j) )
+   enddo
+ enddo
+  
 END_PROVIDER
 
 

src/Hartree_Fock/diagonalize_fock.irp.f

@@ -66,14 +66,14 @@ BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
  do i = 1,elec_alpha_num
   accu = 0.d0
   do j = 1, elec_alpha_num
-   accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j)
+   accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
   diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
  enddo
  do i = elec_alpha_num+1,mo_tot_num
   accu = 0.d0
   do j = 1, elec_alpha_num
-   accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j)
+   accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
   diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
  enddo