minor modifs

mend

plugins/MRPT/MRPT_Utils.main.irp.f

@@ -46,19 +46,6 @@ end
 
 subroutine routine_2
  implicit none
- integer :: i
+ provide electronic_psi_ref_average_value
- do i = 1, n_core_inact_orb
-  print*,fock_core_inactive_total(i,1,1),fock_core_inactive(i)
- enddo
- double precision :: accu
- accu = 0.d0
- do i = 1, n_act_orb
-  integer :: j_act_orb
-  j_act_orb = list_act(i)
-  accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb,1)
-  print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb,1),one_body_dm_mo_beta(j_act_orb,j_act_orb,1) 
- enddo
- print*,'accu = ',accu
-
 end
 

plugins/MRPT_Utils/MRMP2_density.irp.f

@@ -0,0 +1,46 @@
+BEGIN_PROVIDER [double precision, MRMP2_density, (mo_tot_num_align, mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l
+ double precision :: accu, mp2_dm(mo_tot_num)
+ MRMP2_density = one_body_dm_mo
+ call give_2h2p_density(mp2_dm)
+ accu = 0.d0
+ do i = 1, n_virt_orb
+  j = list_virt(i)
+  accu += mp2_dm(j)
+  MRMP2_density(j,j)+= mp2_dm(j)
+ enddo
+
+END_PROVIDER
+
+subroutine give_2h2p_density(mp2_density_diag_alpha_beta)
+ implicit none
+ double precision, intent(out) :: mp2_density_diag_alpha_beta(mo_tot_num)
+ integer :: i,j,k,l,m
+ integer :: iorb,jorb,korb,lorb
+
+ double precision               :: get_mo_bielec_integral
+ double precision               :: direct_int
+ double precision               :: coef_double
+
+ mp2_density_diag_alpha_beta = 0.d0
+   do k = 1, n_virt_orb
+    korb = list_virt(k)
+       do i = 1, n_inact_orb
+        iorb = list_inact(i)
+        do j = 1, n_inact_orb
+         jorb = list_inact(j)
+          do l = 1, n_virt_orb
+              lorb = list_virt(l)
+              direct_int = get_mo_bielec_integral(iorb,jorb,korb,lorb ,mo_integrals_map)
+              coef_double = direct_int/(fock_core_inactive_total_spin_trace(iorb,1) + fock_core_inactive_total_spin_trace(jorb,1) &
+                                     -fock_virt_total_spin_trace(korb,1) - fock_virt_total_spin_trace(lorb,1))
+              mp2_density_diag_alpha_beta(korb) += coef_double * coef_double
+          enddo
+         enddo
+        enddo
+  print*, mp2_density_diag_alpha_beta(korb)
+ enddo
+
+end
+

plugins/MRPT_Utils/mrpt_dress.irp.f

@@ -121,7 +121,8 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
         delta_e(i_state) = 1.d+20
        enddo
       else 
-       call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
+       call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),delta_e)
+!      print*, 'delta_e',delta_e
        !!!!!!!!!!!!! SHIFTED BK 
 !      double precision :: hjj
 !      call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
@@ -129,6 +130,7 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
 !      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       endif
       hij_array(index_i) = hialpha
+!     print*, 'hialpha ',hialpha
       do i_state = 1,N_states
        delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
       enddo

plugins/MRPT_Utils/mrpt_utils.irp.f

@@ -40,7 +40,6 @@
  enddo
  print*, '1h   = ',accu
  
- stop
  ! 1p 
  delta_ij_tmp = 0.d0
  call H_apply_mrpt_1p(delta_ij_tmp,N_det)
@@ -235,7 +234,7 @@ END_PROVIDER
  enddo
  END_PROVIDER 
 
-  BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ]
+  BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_electronic_energy, (N_states_diag) ]
  &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states_diag) ]
  &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ]
   BEGIN_DOC
@@ -279,7 +278,7 @@ END_PROVIDER
      allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
      allocate (eigenvalues(N_det))
      call lapack_diag(eigenvalues,eigenvectors,                      &
-         H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
+         Hmatrix_dressed_pt2_new_symmetrized,size(H_matrix_all_dets,1),N_det)
      CI_electronic_energy(:) = 0.d0
      if (s2_eig) then
        i_state = 0
@@ -288,8 +287,11 @@ END_PROVIDER
        good_state_array = .False.
        call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
          N_det,size(eigenvectors,1))
+!      print*, s2_eigvalues(:)
+       print*,'N_det',N_det
        do j=1,N_det
          ! Select at least n_states states with S^2 values closed to "expected_s2"
+         print*, s2_eigvalues(j),expected_s2
          if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
            i_state +=1
            index_good_state_array(i_state) = j
@@ -303,10 +305,10 @@ END_PROVIDER
          ! Fill the first "i_state" states that have a correct S^2 value
          do j = 1, i_state
            do i=1,N_det
-             CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
+             CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
            enddo
-           CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
+           CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(index_good_state_array(j))
-           CI_eigenvectors_s2(j) = s2_eigvalues(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
@@ -316,10 +318,10 @@ END_PROVIDER
              exit
            endif
            do i=1,N_det
-             CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
+             CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
            enddo
-           CI_electronic_energy(i_state+i_other_state) = eigenvalues(j)
+           CI_dressed_pt2_new_electronic_energy(i_state+i_other_state) = eigenvalues(j)
-           CI_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
+           CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
          enddo
 
        else
@@ -334,10 +336,10 @@ END_PROVIDER
          print*,''
          do j=1,min(N_states_diag,N_det)
            do i=1,N_det
-             CI_eigenvectors(i,j) = eigenvectors(i,j)
+             CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
            enddo
-           CI_electronic_energy(j) = eigenvalues(j)
+           CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j)
-           CI_eigenvectors_s2(j) = s2_eigvalues(j)
+           CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j)
          enddo
        endif
        deallocate(index_good_state_array,good_state_array)
@@ -348,9 +350,9 @@ END_PROVIDER
        ! Select the "N_states_diag" states of lowest energy
        do j=1,min(N_det,N_states_diag)
          do i=1,N_det
-           CI_eigenvectors(i,j) = eigenvectors(i,j)
+           CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
          enddo
-         CI_electronic_energy(j) = eigenvalues(j)
+         CI_dressed_pt2_new_electronic_energy(j) = eigenvalues(j)
        enddo
      endif
      deallocate(eigenvectors,eigenvalues)
@@ -370,7 +372,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
   character*(8)                  :: st
   call write_time(output_determinants)
   do j=1,N_states_diag
-    CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
+    CI_dressed_pt2_new_energy(j) = CI_dressed_pt2_new_electronic_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))

plugins/Psiref_CAS/psi_ref.irp.f

@@ -72,9 +72,19 @@ END_PROVIDER
 &BEGIN_PROVIDER  [double precision, psi_ref_average_value, (N_states)]
  implicit none
  integer :: i,j
- call u_0_H_u_0(electronic_psi_ref_average_value,psi_ref_coef,N_det_ref,psi_ref,N_int,N_states,psi_det_size) 
+ electronic_psi_ref_average_value = psi_energy
  do i = 1, N_states
-  psi_ref_average_value(i) = electronic_psi_ref_average_value(i) + nuclear_repulsion
+  psi_ref_average_value(i) = psi_energy(i) + nuclear_repulsion
  enddo
+ double precision :: accu,hij
+ accu = 0.d0
+ do i = 1, N_det_ref
+  do j = 1, N_det_ref
+   call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
+   accu += psi_ref_coef(i,1) * psi_ref_coef(j,1) * hij
+  enddo
+ enddo
+ electronic_psi_ref_average_value(1) = accu
+ psi_ref_average_value(1) = electronic_psi_ref_average_value(1) + nuclear_repulsion
 
 END_PROVIDER

src/MO_Basis/rotate_mos.irp.f

@@ -0,0 +1,9 @@
+program rotate_mos
+ implicit none
+ integer :: i,j
+ write(*,*)'Which couple of MOs would you like to mix ?'
+ read(5,*)i,j
+ call mix_mo_jk(i,j)
+ call save_mos
+
+end

src/Utils/invert.irp.f

@@ -0,0 +1,19 @@
+subroutine invert_matrix(A,LDA,na,A_inv,LDA_inv)
+implicit none
+double precision, intent(in)    :: A (LDA,na)
+integer, intent(in)             :: LDA, LDA_inv
+integer, intent(in)             :: na
+double precision, intent(out)   :: A_inv (LDA_inv,na)
+
+ double precision :: work(LDA_inv*max(na,64))
+!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: work
+ integer          :: inf
+ integer          :: ipiv(LDA_inv)
+!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
+  integer          :: lwork
+  A_inv(1:na,1:na) = A(1:na,1:na) 
+  call dgetrf(na, na, A_inv, LDA_inv, ipiv, inf )
+  lwork = SIZE(work)
+  call dgetri(na, A_inv, LDA_inv, ipiv, work, lwork, inf )
+end
+