minor modifs in printing

plugins/MRPT_Utils/EZFIO.cfg

@@ -5,3 +5,10 @@ interface: ezfio,provider,ocaml
 default: True
 
 
+[pure_state_specific_mrpt2]
+type: logical
+doc: If true, diagonalize the dressed matrix for each state and do a state following of the initial states
+interface: ezfio,provider,ocaml
+default: True
+
+

plugins/MRPT_Utils/MRPT_Utils.main.irp.f

@@ -10,16 +10,20 @@ end
 
 subroutine routine_3
  implicit none
+ integer :: i
 !provide fock_virt_total_spin_trace
  provide delta_ij 
  
  print *,  'N_det    = ', N_det
  print *,  'N_states = ', N_states
- print *,  'PT2      = ', second_order_pt_new(1)
+ do i = 1, N_States
- print *,  'E        = ', CI_energy(1)
+  print*,'State',i
- print *,  'E+PT2    = ', CI_energy(1)+second_order_pt_new(1)
+  write(*,'(A12,X,I3,A3,XX,F16.10)')    '  PT2       ', i,' = ', second_order_pt_new(i)
- print *,'****** DIAGONALIZATION OF DRESSED MATRIX ******'
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E         ', i,' = ', CI_energy(i)
- print *,  'E dressed= ', CI_dressed_pt2_new_energy(1)
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E+PT2     ', i,' = ', CI_energy(i)+second_order_pt_new(i)
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  S^2       ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
+ enddo 
 
 end
 

plugins/MRPT_Utils/energies_cas.irp.f

@@ -241,13 +241,13 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
                norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
        call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
                norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states)
-       if(orb_i == orb_j .and. ispin .ne. jspin)then  
+      !if(orb_i == orb_j .and. ispin .ne. jspin)then  
         call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
         one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target)  -   energies(state_target)
-       else
+      !else
-        call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
+      ! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
-        one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target)  -   energies(state_target)
+      ! one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target)  -   energies(state_target)
-       endif
+      !endif
       enddo
     enddo
    enddo
@@ -527,7 +527,7 @@ END_PROVIDER
 
  double precision :: thresh_norm
   
- thresh_norm = 1.d-10
+ thresh_norm = 1.d-20
 
 
 
@@ -552,14 +552,7 @@ END_PROVIDER
        call debug_det(psi_in_out,N_int)
        print*, 'pb, i_ok ne 0 !!!'
       endif
-!     call i_H_j_no_k_operators_from_act(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij_test)
       call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
-!     if(i==1.and.dabs(hij)>1.d-8)then
-!     if(dabs(hij)>1.d-8)then
-!      print*, ispin,vorb,iorb
-!      print*, i,hij,hij_test
-!      pause
-!     endif
       do j = 1, n_states
         double precision ::  coef,contrib
         coef = psi_coef(i,j) !* psi_coef(i,j)
@@ -635,7 +628,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
 
  double precision :: thresh_norm
   
- thresh_norm = 1.d-10
+ thresh_norm = 1.d-20
 
  do aorb = 1,n_act_orb
   orb_a = list_act(aorb)

plugins/MRPT_Utils/mrpt_dress.irp.f

@@ -84,7 +84,11 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
       do i_state = 1, N_states
        coef_array(i_state) = psi_coef(index_i,i_state)
       enddo
-      call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
+      if(dabs(hialpha).le.1.d-10)then
+       delta_e = 1.d+20
+      else
+       call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
+      endif
       hij_array(index_i) = hialpha
       call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
 !     phase_array(index_i) = phase

plugins/MRPT_Utils/mrpt_utils.irp.f

@@ -121,8 +121,8 @@
 
  ! 1h2p   
  delta_ij_tmp = 0.d0
-!call give_1h2p_contrib(delta_ij_tmp)
+ call give_1h2p_contrib(delta_ij_tmp)
- call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
+!!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
  accu = 0.d0
  do i_state = 1, N_states
  do i = 1, N_det
@@ -137,8 +137,8 @@
 
  ! 2h1p   
  delta_ij_tmp = 0.d0
-!call give_2h1p_contrib(delta_ij_tmp)
+ call give_2h1p_contrib(delta_ij_tmp)
- call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
+!!!! call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
  accu = 0.d0
  do i_state = 1, N_states
  do i = 1, N_det
@@ -152,19 +152,19 @@
  print*, '2h1p = ',accu
 
  ! 2h2p   
-!delta_ij_tmp = 0.d0
+ delta_ij_tmp = 0.d0
-!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det)
+!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det)
-!accu = 0.d0
+ accu = 0.d0
-!do i_state = 1, N_states
+ do i_state = 1, N_states
-!do i = 1, N_det
+ do i = 1, N_det
-! do j = 1, N_det
+  do j = 1, N_det
-!  accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
+   accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
-!  delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
-! enddo
+  enddo
-!enddo
+ enddo
-!second_order_pt_new_2h2p(i_state) = accu(i_state) 
+ second_order_pt_new_2h2p(i_state) = accu(i_state) 
-!enddo
+ enddo
-!print*, '2h2p = ',accu
+ print*, '2h2p = ',accu
 
  double precision :: contrib_2h2p(N_states)
  call give_2h2p(contrib_2h2p)
@@ -236,13 +236,15 @@ END_PROVIDER
   logical, allocatable  :: good_state_array(:)
   double precision, allocatable :: s2_values_tmp(:)
   integer :: i_other_state
-  double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+  double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:), hmatrix_tmp(:,:)
   integer :: i_state
   double precision :: s2,e_0
   integer                        :: i,j,k
   double precision, allocatable :: s2_eigvalues(:)
   double precision, allocatable :: e_array(:)
   integer, allocatable :: iorder(:)
+  double precision :: overlap(N_det)
+  double precision, allocatable :: psi_tmp(:)
   
   ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors 
   do j=1,min(N_states,N_det)
@@ -265,82 +267,131 @@ END_PROVIDER
   else if (diag_algorithm == "Lapack") then
      allocate (eigenvectors(N_det,N_det))
      allocate (eigenvalues(N_det))
-     call lapack_diag(eigenvalues,eigenvectors,                      &
+     if(pure_state_specific_mrpt2)then
-         Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det)
+      allocate (hmatrix_tmp(N_det,N_det))
-     CI_electronic_dressed_pt2_new_energy(:) = 0.d0
+      allocate (iorder(N_det))
-     if (s2_eig) then
+      allocate (psi_tmp(N_det))
-       i_state = 0
+       print*,''
-       allocate (s2_eigvalues(N_det))
+       print*,'***************************'
-       allocate(index_good_state_array(N_det),good_state_array(N_det))
+      do i_state = 1, N_states !! Big loop over states
-       good_state_array = .False.
+       print*,''
-       call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
+       print*,'Diagonalizing with the dressing for state',i_state
-         N_det,size(eigenvectors,1))
+       do i = 1, N_det
-       do j=1,N_det
+        do j = 1, N_det
-         ! Select at least n_states states with S^2 values closed to "expected_s2"
+         hmatrix_tmp(j,i) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state)
-         if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
+        enddo
-           i_state += 1
-           index_good_state_array(i_state) = j
-           good_state_array(j) = .True.
-         endif
-         if (i_state==N_states) then
-           exit
-         endif
        enddo
-       if (i_state /= 0) then
+       call lapack_diag(eigenvalues,eigenvectors,                      &
-         ! Fill the first "i_state" states that have a correct S^2 value
+           Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det)
-         do j = 1, i_state
+       write(*,'(A86)')'Looking for the most overlapping state within all eigenvectors of the dressed matrix'
-           do i=1,N_det
+       print*,''
-             CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
+       print*,'Calculating the overlap for ...'
-           enddo
+       do i = 1, N_det
-           CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j))
+        overlap(i) = 0.d0
-           CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
+        iorder(i) = i
-         enddo
+        print*,'eigenvector',i
-         i_other_state = 0
+        do j = 1, N_det
-         do j = 1, N_det
+         overlap(i)+= psi_coef(j,i_state) * eigenvectors(j,i)
-           if(good_state_array(j))cycle
+        enddo
-           i_other_state +=1
+        overlap(i) = -dabs(overlap(i))
-           if(i_state+i_other_state.gt.n_states_diag)then
+        print*,'energy  = ',eigenvalues(i) + nuclear_repulsion
-             exit
+        print*,'overlap = ',dabs(overlap(i))
-           endif
-           do i=1,N_det
-             CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
-           enddo
-           CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j)
-           CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
-         enddo
-         
-       else
-         print*,''
-         print*,'!!!!!!!!   WARNING  !!!!!!!!!'
-         print*,'  Within the ',N_det,'determinants selected'
-         print*,'  and the ',N_states_diag,'states requested'
-         print*,'  We did not find any state with S^2 values close to ',expected_s2
-         print*,'  We will then set the first N_states eigenvectors of the H matrix'
-         print*,'  as the CI_dressed_pt2_new_eigenvectors'
-         print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
-         print*,''
-         do j=1,min(N_states_diag,N_det)
-           do i=1,N_det
-             CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
-           enddo
-           CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
-           CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j)
-         enddo
-       endif
-       deallocate(index_good_state_array,good_state_array)
-       deallocate(s2_eigvalues)
-     else
-       call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,&
-          min(N_det,N_states_diag),size(eigenvectors,1))
-       ! Select the "N_states_diag" states of lowest energy
-       do j=1,min(N_det,N_states_diag)
-         do i=1,N_det
-           CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
-         enddo
-         CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
        enddo
+       print*,''
+       print*,'Sorting the eigenvectors per overlap'
+       call dsort(overlap,iorder,n_states)
+       print*,''
+       print*,'The most overlapping state is the ',iorder(1)
+       print*,'with the overlap of ',dabs(overlap(1))
+       print*,'and an energy of    ',eigenvalues(iorder(1)) + nuclear_repulsion
+       print*,'Calculating the S^2 value ...'
+       do i=1,N_det
+         CI_dressed_pt2_new_eigenvectors(i,i_state) = eigenvectors(i,iorder(1))
+         psi_tmp(i) = eigenvectors(i,iorder(1))
+       enddo
+       CI_electronic_dressed_pt2_new_energy(i_state) = eigenvalues(iorder(1))
+       call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det,psi_det,N_int,1,N_det)
+       print*,'S^2      = ', CI_dressed_pt2_new_eigenvectors_s2(i_state)
+      enddo
+     else 
+      call lapack_diag(eigenvalues,eigenvectors,                      &
+          Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det)
+      CI_electronic_dressed_pt2_new_energy(:) = 0.d0
+      if (s2_eig) then
+        i_state = 0
+        allocate (s2_eigvalues(N_det))
+        allocate(index_good_state_array(N_det),good_state_array(N_det))
+        good_state_array = .False.
+        call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
+          N_det,size(eigenvectors,1))
+        do j=1,N_det
+          ! Select at least n_states states with S^2 values closed to "expected_s2"
+          print*, eigenvalues(j)+nuclear_repulsion, s2_eigvalues(j)
+          if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
+            i_state += 1
+            index_good_state_array(i_state) = j
+            good_state_array(j) = .True.
+          endif
+          if (i_state==N_states) then
+            exit
+          endif
+        enddo
+        if (i_state /= 0) then
+          ! Fill the first "i_state" states that have a correct S^2 value
+          do j = 1, i_state
+            do i=1,N_det
+              CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
+            enddo
+            CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j))
+            CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
+          enddo
+          i_other_state = 0
+          do j = 1, N_det
+            if(good_state_array(j))cycle
+            i_other_state +=1
+            if(i_state+i_other_state.gt.n_states)then
+              exit
+            endif
+            do i=1,N_det
+              CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
+            enddo
+            CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j)
+            CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
+          enddo
+          
+        else
+          print*,''
+          print*,'!!!!!!!!   WARNING  !!!!!!!!!'
+          print*,'  Within the ',N_det,'determinants selected'
+          print*,'  and the ',N_states_diag,'states requested'
+          print*,'  We did not find any state with S^2 values close to ',expected_s2
+          print*,'  We will then set the first N_states eigenvectors of the H matrix'
+          print*,'  as the CI_dressed_pt2_new_eigenvectors'
+          print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
+          print*,''
+          do j=1,min(N_states_diag,N_det)
+            do i=1,N_det
+              CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
+            enddo
+            CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
+            CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j)
+          enddo
+        endif
+        deallocate(index_good_state_array,good_state_array)
+        deallocate(s2_eigvalues)
+      else
+        call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,&
+           min(N_det,N_states_diag),size(eigenvectors,1))
+        ! Select the "N_states_diag" states of lowest energy
+        do j=1,min(N_det,N_states)
+          do i=1,N_det
+            CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
+          enddo
+          CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j)
+        enddo
+      endif
+      deallocate(eigenvectors,eigenvalues)
      endif
-     deallocate(eigenvectors,eigenvalues)
 
   endif
 
@@ -361,7 +412,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
     CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
     write(st,'(I4)') j
     call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))
-    call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
+    call write_double(output_determinants, CI_dressed_pt2_new_eigenvectors_s2(j) ,'S^2 of state '//trim(st))
   enddo
 
 END_PROVIDER

plugins/MRPT_Utils/print_1h2p.irp.f

@@ -8,8 +8,10 @@ end
 subroutine routine_2
  implicit none
  integer :: i,j
- do i =1, n_inact_orb
+ do i =1, n_act_orb
-  write(*,'(100(F16.10,X))')one_anhil_one_creat_inact_virt(i,:,1)
+!do i =1, 2
+  write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(i,:,:,:,1)
+! write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(1,4,1,2,1)
  enddo
 
 

plugins/MRPT_Utils/special_hij.irp.f

@@ -1,183 +0,0 @@
-
-
-subroutine i_H_j_no_k_operators_from_act(key_i,key_j,Nint,hij)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Returns <i|H|j> where i and j are determinants
-  END_DOC
-  integer, intent(in)            :: Nint
-  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
-  double precision, intent(out)  :: hij
-
-  integer                        :: exc(0:2,2,2)
-  integer                        :: degree
-  double precision               :: get_mo_bielec_integral, phase
-  integer                        :: m,n,p,q
-  integer                        :: i,j,k
-  integer                        :: occ(Nint*bit_kind_size,2)
-  double precision               :: diag_H_mat_elem
-  integer                        :: n_occ_ab(2)
-  logical                        :: has_mipi(Nint*bit_kind_size)
-  double precision               :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size), miip_other(Nint*bit_kind_size)
-  PROVIDE mo_bielec_integrals_in_map mo_integrals_map
-
-  ASSERT (Nint > 0)
-  ASSERT (Nint == N_int)
-  ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
-  ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
-  ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
-  ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
-
-  hij = 0.d0
-  !DIR$ FORCEINLINE
-  call get_excitation_degree(key_i,key_j,degree,Nint)
-  select case (degree)
-    case (2)
-      call get_double_excitation(key_i,key_j,exc,phase,Nint)
-      if (exc(0,1,1) == 1) then
-        ! Mono alpha, mono beta
-        hij = phase*get_mo_bielec_integral(                          &
-            exc(1,1,1),                                              &
-            exc(1,1,2),                                              &
-            exc(1,2,1),                                              &
-            exc(1,2,2) ,mo_integrals_map)
-      else if (exc(0,1,1) == 2) then
-        ! Double alpha
-        hij = phase*(get_mo_bielec_integral(                         &
-            exc(1,1,1),                                              &
-            exc(2,1,1),                                              &
-            exc(1,2,1),                                              &
-            exc(2,2,1) ,mo_integrals_map) -                          &
-            get_mo_bielec_integral(                                  &
-            exc(1,1,1),                                              &
-            exc(2,1,1),                                              &
-            exc(2,2,1),                                              &
-            exc(1,2,1) ,mo_integrals_map) )
-      else if (exc(0,1,2) == 2) then
-        ! Double beta
-        hij = phase*(get_mo_bielec_integral(                         &
-            exc(1,1,2),                                              &
-            exc(2,1,2),                                              &
-            exc(1,2,2),                                              &
-            exc(2,2,2) ,mo_integrals_map) -                          &
-            get_mo_bielec_integral(                                  &
-            exc(1,1,2),                                              &
-            exc(2,1,2),                                              &
-            exc(2,2,2),                                              &
-            exc(1,2,2) ,mo_integrals_map) )
-      endif
-    case (1)
-      call get_mono_excitation(key_i,key_j,exc,phase,Nint)
-      !DIR$ FORCEINLINE
-      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
-      has_mipi = .False.
-      logical :: is_i_in_active
-      double precision :: accu_a, accu_b, accu_core
-      accu_a = 0.d0
-      accu_b = 0.d0
-      if (exc(0,1,1) == 1) then
-        ! Mono alpha
-        m = exc(1,1,1)
-        p = exc(1,2,1)
-        do k = 1, elec_alpha_num
-          i = occ(k,1)
-          if (.not.has_mipi(i)) then
-            mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
-            has_mipi(i) = .True.
-          endif
-          if(.not.is_i_in_active(i))then
-           miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-          else 
-!          print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-           miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-           accu_a += miip(i)
-          endif
-        enddo
-        do k = 1, elec_beta_num
-          i = occ(k,2)
-          if (.not.has_mipi(i)) then
-            mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
-            has_mipi(i) = .True.
-          endif
-          if(.not.is_i_in_active(i))then
-           miip_other(i) = 0.d0
-          else 
-!          print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-           miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-           accu_b += miip(i)
-          endif
-        enddo
-!       print*, accu_a,accu_b,accu_a + accu_b
-        accu_a = 0.d0
-        accu_b = 0.d0
-        accu_core = mo_mono_elec_integral(m,p)
-
-        do k = 1, elec_alpha_num
-          hij = hij + mipi(occ(k,1))  - miip(occ(k,1))
-          accu_Core += mipi(occ(k,1))
-          if(is_i_in_active(occ(k,1)))then
-           accu_a += miip(occ(k,1))
-          else 
-           accu_Core -= miip(occ(k,1))
-          endif
-        enddo
-!       print*, hij,accu_core
-        do k = 1, elec_beta_num
-          hij = hij + mipi(occ(k,2))  - miip_other(occ(k,2))
-          accu_Core += mipi(occ(k,2))
-          if(is_i_in_active(occ(k,2)))then
-           accu_b += miip_other(occ(k,2))
-          else 
-           accu_Core -= miip_other(occ(k,2))
-          endif
-        enddo
-!       print*, hij,accu_core,accu_core - accu_a - accu_b
-!       print*, accu_a,accu_b,accu_a + accu_b
-!       pause
-
-      else
-        ! Mono beta
-        m = exc(1,1,2)
-        p = exc(1,2,2)
-        do k = 1, elec_beta_num
-          i = occ(k,2)
-          if (.not.has_mipi(i)) then
-            mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
-            has_mipi(i) = .True.
-          endif
-          if(.not.is_i_in_active(i))then
-           miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-          else 
-           miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-          endif
-        enddo
-        do k = 1, elec_alpha_num
-          i = occ(k,1)
-          if (.not.has_mipi(i)) then
-            mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
-            has_mipi(i) = .True.
-          endif
-          if(.not.is_i_in_active(i))then
-           miip_other(i) = 0.d0
-          else 
-           miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
-          endif
-        enddo
-        do k = 1, elec_alpha_num
-          hij = hij + mipi(occ(k,1)) - miip_other(occ(k,1))
-        enddo
-        do k = 1, elec_beta_num
-          hij = hij + mipi(occ(k,2)) - miip(occ(k,2))
-        enddo
-
-      endif
-      hij = phase*(hij + mo_mono_elec_integral(m,p))
-
-    case (0)
-      hij = diag_H_mat_elem(key_i,Nint)
-  end select
-end
-
-
-

plugins/loc_cele/loc_exchange_int.irp.f

@@ -13,21 +13,21 @@ program loc_int
   iorb = list_core_inact(i)
   exchange_int = 0.d0
   iorder = 0
-  print*,''
   if(list_core_inact_check(iorb) == .False.)cycle
   do j = i+1, n_core_inact_orb
    jorb = list_core_inact(j)
    iorder(jorb) = jorb
-   exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   if(list_core_inact_check(jorb) == .False.)then
+    exchange_int(jorb) = 0.d0
+   else 
+    exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   endif
   enddo
   n_rot += 1
   call dsort(exchange_int,iorder,mo_tot_num)
   indices(n_rot,1) = iorb
   indices(n_rot,2) = iorder(1)
   list_core_inact_check(iorder(1)) = .False.
-  print*,indices(n_rot,1),indices(n_rot,2)
-  print*,''
-  print*,''
  enddo
   print*,'****************************'
   print*,'-+++++++++++++++++++++++++'
@@ -45,21 +45,21 @@ program loc_int
   iorb = list_act(i)
   exchange_int = 0.d0
   iorder = 0
-  print*,''
   if(list_core_inact_check(iorb) == .False.)cycle
   do j = i+1, n_act_orb
    jorb = list_act(j)
    iorder(jorb) = jorb
-   exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   if(list_core_inact_check(jorb) == .False.)then
+    exchange_int(jorb) = 0.d0
+   else 
+    exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   endif
   enddo
   n_rot += 1
   call dsort(exchange_int,iorder,mo_tot_num)
   indices(n_rot,1) = iorb
   indices(n_rot,2) = iorder(1)
   list_core_inact_check(iorder(1)) = .False.
-  print*,indices(n_rot,1),indices(n_rot,2)
-  print*,''
-  print*,''
  enddo
   print*,'****************************'
   print*,'-+++++++++++++++++++++++++'
@@ -82,16 +82,17 @@ program loc_int
   do j = i+1, n_virt_orb
    jorb = list_virt(j)
    iorder(jorb) = jorb
-   exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   if(list_core_inact_check(jorb) == .False.)then
+    exchange_int(jorb) = 0.d0
+   else 
+    exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   endif
   enddo
   n_rot += 1
   call dsort(exchange_int,iorder,mo_tot_num)
   indices(n_rot,1) = iorb
   indices(n_rot,2) = iorder(1)
   list_core_inact_check(iorder(1)) = .False.
-  print*,indices(n_rot,1),indices(n_rot,2)
-  print*,''
-  print*,''
  enddo
   print*,'****************************'
   print*,'-+++++++++++++++++++++++++'

src/Determinants/slater_rules.irp.f

@@ -2168,9 +2168,27 @@ subroutine H_u_0_stored(v_0,u_0,hmatrix,sze)
   double precision, intent(in)   :: u_0(sze)
   v_0 = 0.d0
   call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
-
 end
 
+subroutine H_s2_u_0_stored(v_0,u_0,hmatrix,s2matrix,sze)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes v_0 = H|u_0>
+  !
+  ! n : number of determinants
+  !
+  ! uses the big_matrix_stored array
+  END_DOC
+  integer, intent(in)            :: sze
+  double precision, intent(in)   :: hmatrix(sze,sze),s2matrix(sze,sze)
+  double precision, intent(out)  :: v_0(sze)
+  double precision, intent(in)   :: u_0(sze)
+  v_0 = 0.d0
+  call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
+end
+
+
 subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze)
   use bitmasks
   implicit none