added pt2_tc and pt2_tc_cisd

external/qp2-dependencies

@@ -1 +1 @@
-Subproject commit 6e23ebac001acae91d1c762ca934e09a9b7d614a
+Subproject commit e0d0e02e9f5ece138d1520106954a881ab0b8db2

src/cipsi_tc_bi_ortho/pt2.irp.f

@@ -1,4 +1,4 @@
-subroutine pt2_tc_bi_ortho
+subroutine tc_pt2
   use selection_types
   implicit none
   BEGIN_DOC
@@ -15,7 +15,7 @@ subroutine pt2_tc_bi_ortho
   double precision, external :: memory_of_double
   double precision :: correlation_energy_ratio,E_denom,E_tc,norm
   double precision, allocatable :: ept2(:), pt1(:),extrap_energy(:)
-  PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
+  PROVIDE H_apply_buffer_allocated distributed_davidson 
 
   print*,'Diagonal elements of the Fock matrix '
   do i = 1, mo_num
@@ -44,24 +44,14 @@ subroutine pt2_tc_bi_ortho
   pt2_data % overlap= 0.d0
   pt2_data % variance = huge(1.e0)
 
-  if (s2_eig) then
+  !!!! WARNING  !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION
-    call make_s2_eigenfunction
+!  if (s2_eig) then
-  endif
+!    call make_s2_eigenfunction
+!  endif
   print_pt2 = .False.
   call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
 !  call routine_save_right
 
-  if (N_det > N_det_max) then
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det)
-    psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states)
-    N_det = N_det_max
-    soft_touch N_det psi_det psi_coef
-    if (s2_eig) then
-      call make_s2_eigenfunction
-    endif
-    print_pt2 = .False.
-    call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
-  endif
 
   allocate(ept2(1000),pt1(1000),extrap_energy(100))
 
@@ -71,18 +61,11 @@ subroutine pt2_tc_bi_ortho
 ! soft_touch thresh_it_dav
 
   print_pt2 = .True.
-  to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
-  to_select = max(N_states_diag, to_select)
-
-  E_denom = E_tc ! TC Energy of the current wave function 
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
   call pt2_alloc(pt2_data, N_states)
   call pt2_alloc(pt2_data_err, N_states)
-  call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
+  call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection
-
-  N_iter += 1
-
   call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
 
 end

src/fci_tc_bi/pt2_tc.irp.f

@@ -0,0 +1,31 @@
+program tc_pt2_prog
+  implicit none
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid 
+  pruning = -1.d0
+  touch pruning
+!  pt2_relative_error = 0.01d0
+!  touch pt2_relative_error
+  call run_pt2_tc
+
+end
+
+
+subroutine run_pt2_tc
+
+ implicit none
+
+  PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e
+  if(elec_alpha_num+elec_beta_num.ge.3)then
+    if(three_body_h_tc)then
+      call provide_all_three_ints_bi_ortho
+    endif
+  endif
+  ! ---
+
+ call tc_pt2
+
+
+end

src/kohn_sham/print_mos.irp.f

@@ -3,7 +3,7 @@ program print_mos
  integer :: i,nx
  double precision :: r(3), xmax, dx, accu
  double precision, allocatable :: mos_array(:)
- double precision:: alpha,envelop
+ double precision:: alpha,envelop,dm_a,dm_b
  allocate(mos_array(mo_num))
  xmax = 5.d0
  nx = 1000
@@ -11,20 +11,14 @@ program print_mos
  r = 0.d0
  alpha = 0.5d0
  do i = 1, nx
+  call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
   call give_all_mos_at_r(r,mos_array)
   accu = mos_array(3)**2+mos_array(4)**2+mos_array(5)**2
   accu = dsqrt(accu)
   envelop = (1.d0 - dexp(-alpha * r(3)**2))
-  write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, envelop
+  write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, dm_a+dm_b, envelop
   r(3) += dx
  enddo
 
 end
 
-double precision function f_mu(x)
- implicit none
- double precision, intent(in) :: x
- 
-
-
-end

src/tc_bi_ortho/e_corr_bi_ortho.irp.f

@@ -45,6 +45,9 @@
 &BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj ]
 &BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth ]
 &BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth ]
+&BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj_abs ]
+&BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth_abs ]
+&BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth_abs ]
  implicit none 
  integer :: i,degree
  double precision :: hmono,htwoe,hthree,htilde_ij
@@ -57,13 +60,15 @@
   call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
   if(degree == 1)then
    e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
+   e_corr_single_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1))
   else if(degree == 2)then
    e_corr_double_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
-!   print*,'coef_wf , e_cor',reigvec_tc_bi_orth(i,1)/reigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
+   e_corr_double_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1))
   endif
  enddo
  e_corr_bi_orth_proj = e_corr_single_bi_orth + e_corr_double_bi_orth
  e_corr_bi_orth = eigval_right_tc_bi_orth(1) - e_tilde_bi_orth_00
+ e_corr_bi_orth_proj_abs = e_corr_single_bi_orth_abs + e_corr_double_bi_orth_abs
  END_PROVIDER 
 
  BEGIN_PROVIDER [ double precision, e_tc_left_right ]

src/tc_bi_ortho/pt2_tc_cisd.irp.f

@@ -0,0 +1,43 @@
+program pt2_tc_cisd
+
+  BEGIN_DOC
+  !
+  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together 
+  !        with the energy. Saves the left-right wave functions at the end. 
+  !
+  END_DOC
+
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  read_wf = .True.
+  touch read_wf
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  print*, ' nb of states = ', N_states
+  print*, ' nb of det    = ', N_det
+
+  call routine
+end
+
+subroutine routine
+ implicit none
+ integer :: i,h1,p1,h2,p2,s1,s2
+ double precision :: h0i,hi0,e00,ei,delta_e
+ double precision :: norm,e_corr,coef
+ norm = 0.d0 
+ e_corr = 0.d0
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,1), N_int, e00) 
+ do i = 2, N_det
+  call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0) 
+  call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i) 
+  call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, ei) 
+  delta_e = e00 - ei
+  coef = hi0/delta_e
+  norm += coef*coef
+  e_corr += dabs(coef* h0i)
+ enddo
+ print*,'e_corr = ',e_corr
+ print*,'norm   = ',norm
+
+end

src/tc_bi_ortho/slater_tc_opt_diag.irp.f

@@ -156,7 +156,7 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
     htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
   enddo
 
-  if(three_body_h_tc)then
+  if(three_body_h_tc.and.elec_num.gt.2)then
    !!!!! 3-e part 
    !! same-spin/same-spin
    do j = 1, na
@@ -243,7 +243,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
     htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
   enddo
 
-  if(three_body_h_tc)then
+  if(three_body_h_tc.and.elec_num.gt.2)then
    !!!!! 3-e part 
    !! same-spin/same-spin
    do j = 1, na

src/tc_bi_ortho/slater_tc_opt_double.irp.f

@@ -41,15 +41,15 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
   if(s1.ne.s2)then
    ! opposite spin two-body 
     htwoe  = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) 
-    if(three_body_h_tc)then
+    if(three_body_h_tc.and.elec_num.gt.2)then
      if(.not.double_normal_ord)then
       if(degree_i>degree_j)then
        call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
       else
        call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
       endif
-     elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then
+     elseif(double_normal_ord.and.elec_num.gt.2)then
-      htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ???
+      htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)
      endif
     endif
   else
@@ -58,16 +58,16 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
    htwoe  = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)  
    ! exchange terms 
    htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) 
-   if(three_body_h_tc)then
+   if(three_body_h_tc.and.elec_num.gt.2)then
     if(.not.double_normal_ord)then
      if(degree_i>degree_j)then
       call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
      else
       call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
      endif
-    elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then
+    elseif(double_normal_ord.and.elec_num.gt.2)then
-      htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ???
+      htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)
-      htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ???
+      htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)
     endif
    endif
   endif

src/tc_bi_ortho/slater_tc_opt_single.irp.f

@@ -106,7 +106,7 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
   htwoe -= buffer_x(i)
  enddo
  hthree = 0.d0
- if (three_body_h_tc)then
+ if (three_body_h_tc.and.elec_num.gt.2)then
   call three_comp_fock_elem(key_i,h,p,spin,hthree)
  endif
 

src/tc_bi_ortho/tc_bi_ortho.irp.f

@@ -57,8 +57,11 @@ subroutine routine_diag()
     print*,'***'                      
     print*,'e_corr_bi_orth            = ',e_corr_bi_orth
     print*,'e_corr_bi_orth_proj       = ',e_corr_bi_orth_proj
+    print*,'e_corr_bi_orth_proj_abs   = ',e_corr_bi_orth_proj_abs
     print*,'e_corr_single_bi_orth     = ',e_corr_single_bi_orth
     print*,'e_corr_double_bi_orth     = ',e_corr_double_bi_orth
+    print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs
+    print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs
     print*,'Left/right eigenvectors'
     do i = 1,N_det
       write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)