added some stuffs for TC-CASSCF

external/ezfio

@@ -1 +1 @@
-Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
+Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c

external/irpf90

@@ -1 +1 @@
-Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102
+Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271

src/casscf_tc_bi/det_manip.irp.f

@@ -0,0 +1,125 @@
+use bitmasks
+
+subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
+      ispin,phase,ierr)
+  BEGIN_DOC
+  ! we create the mono-excitation, and determine, if possible,
+  ! the phase and the number in the list of determinants
+  END_DOC
+  implicit none
+  integer(bit_kind)              :: key1(N_int,2),key2(N_int,2)
+  integer(bit_kind), allocatable :: keytmp(:,:)
+  integer                        :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin
+  real*8                         :: phase
+  logical                        :: found
+  allocate(keytmp(N_int,2))
+  
+  nu=-1
+  phase=1.D0
+  ierr=0
+  call det_copy(key1,key2,N_int)
+  !        write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin
+  !        call print_det(key2,N_int)
+  call do_single_excitation(key2,ihole,ipart,ispin,ierr)
+  !        write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin
+  !        call print_det(key2,N_int)
+  !        write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr
+  if (ierr.eq.1) then
+    ! excitation is possible
+    ! get the phase
+    call get_single_excitation(key1,key2,exc,phase,N_int)
+    ! get the number in the list
+    found=.false.
+    nu=0
+
+    !TODO BOTTLENECK
+    do while (.not.found)
+      nu+=1
+      if (nu.gt.N_det) then
+        ! the determinant is possible, but not in the list
+        found=.true.
+        nu=-1
+      else
+        call det_extract(keytmp,nu,N_int)
+        integer                        :: i,ii
+        found=.true.
+        do ii=1,2
+          do i=1,N_int
+            if (keytmp(i,ii).ne.key2(i,ii)) then
+              found=.false.
+            end if
+          end do
+        end do
+      end if
+    end do
+  end if
+  !
+  ! we found the new string, the phase, and possibly the number in the list
+  !
+end subroutine do_signed_mono_excitation
+
+subroutine det_extract(key,nu,Nint)
+  BEGIN_DOC
+  ! extract a determinant from the list of determinants
+  END_DOC
+  implicit none
+  integer                        :: ispin,i,nu,Nint
+  integer(bit_kind)              :: key(Nint,2)
+  do ispin=1,2
+    do i=1,Nint
+      key(i,ispin)=psi_det(i,ispin,nu)
+    end do
+  end do
+end subroutine det_extract
+
+subroutine det_copy(key1,key2,Nint)
+  use bitmasks ! you need to include the bitmasks_module.f90 features
+  BEGIN_DOC
+  ! copy a determinant from key1 to key2
+  END_DOC
+  implicit none
+  integer                        :: ispin,i,Nint
+  integer(bit_kind)              :: key1(Nint,2),key2(Nint,2)
+  do ispin=1,2
+    do i=1,Nint
+      key2(i,ispin)=key1(i,ispin)
+    end do
+  end do
+end subroutine det_copy
+
+subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2      &
+      ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr)
+  BEGIN_DOC
+  ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q)
+  ! we may create two determinants as result
+  !
+  END_DOC
+  implicit none
+  integer(bit_kind)              :: key_in(N_int,2),key_out1(N_int,2)
+  integer(bit_kind)              :: key_out2(N_int,2)
+  integer                        :: ihole,ipart,ierr,jerr,nu1,nu2
+  integer                        :: ispin
+  real*8                         :: phase1,phase2
+  
+  !        write(6,*) ' applying E_',ipart,ihole,' on determinant '
+  !        call print_det(key_in,N_int)
+  
+  ! spin alpha
+  ispin=1
+  call do_signed_mono_excitation(key_in,key_out1,nu1,ihole           &
+      ,ipart,ispin,phase1,ierr)
+  !        if (ierr.eq.1) then
+  !         write(6,*) ' 1 result is ',nu1,phase1
+  !         call print_det(key_out1,N_int)
+  !        end if
+  ! spin beta
+  ispin=2
+  call do_signed_mono_excitation(key_in,key_out2,nu2,ihole           &
+      ,ipart,ispin,phase2,jerr)
+  !        if (jerr.eq.1) then
+  !         write(6,*) ' 2 result is ',nu2,phase2
+  !         call print_det(key_out2,N_int)
+  !        end if
+  
+end subroutine do_spinfree_mono_excitation
+

src/casscf_tc_bi/grad_dm.irp.f

@@ -193,7 +193,7 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
    r = list_act(rr)
    res_r_inact_test += -tc_transition_matrix_mo(r,t,1,1) * & 
    (2.d0 * mo_bi_ortho_tc_two_e(r,j,a,j) - mo_bi_ortho_tc_two_e(r,j,j,a))
-   res_l_inact_test += -tc_transition_matrix_mo(t,r,1,1) * & 
+   res_l_inact_test -= -tc_transition_matrix_mo(t,r,1,1) * & 
    (2.d0 * mo_bi_ortho_tc_two_e(a,j,r,j) - mo_bi_ortho_tc_two_e(j,a,r,j))
   enddo
  enddo
@@ -207,7 +207,7 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
     u = list_act(uu)
     res_r_act_test += - (mo_bi_ortho_tc_two_e(v,r,u,a) * tc_two_rdm(r,v,t,u) & 
                         +mo_bi_ortho_tc_two_e(v,r,a,u) * tc_two_rdm(r,v,u,t))
-    res_l_act_test += - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) & 
+    res_l_act_test -= - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) & 
                         +mo_bi_ortho_tc_two_e(a,u,v,r) * tc_two_rdm(u,t,r,v))
    enddo
   enddo

src/casscf_tc_bi/test_tc_casscf.irp.f

@@ -0,0 +1,252 @@
+program tc_bi_ortho
+
+  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_i_h_psi
+!  call routine_grad
+ call routine_grad_num_dm_one_body
+end
+
+subroutine routine_i_h_psi
+ implicit none
+ integer :: i,j
+ double precision :: i_H_chi_array(0:3,N_states),i_H_phi_array(0:3,N_states)
+ double precision               :: hmono, htwoe, hthree, htot
+ double precision               :: accu_l_hmono, accu_l_htwoe, accu_l_hthree, accu_l_htot
+ double precision               :: accu_r_hmono, accu_r_htwoe, accu_r_hthree, accu_r_htot
+ double precision               :: test_l_hmono, test_l_htwoe, test_l_hthree, test_l_htot
+ double precision               :: test_r_hmono, test_r_htwoe, test_r_hthree, test_r_htot
+
+ test_l_hmono = 0.d0
+ test_l_htwoe = 0.d0
+ test_l_hthree= 0.d0
+ test_l_htot  = 0.d0
+ test_r_hmono = 0.d0
+ test_r_htwoe = 0.d0
+ test_r_hthree= 0.d0
+ test_r_htot  = 0.d0
+
+ do i = 1, N_det
+   call i_H_tc_psi_phi(psi_det(1,1,i),psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,&
+                       N_int,N_det,N_det,N_states,i_H_chi_array,i_H_phi_array)
+   accu_l_hmono = 0.d0
+   accu_l_htwoe = 0.d0
+   accu_l_hthree= 0.d0
+   accu_l_htot  = 0.d0
+   accu_r_hmono = 0.d0
+   accu_r_htwoe = 0.d0
+   accu_r_hthree= 0.d0
+   accu_r_htot  = 0.d0
+   do j = 1, N_det
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
+    accu_l_hmono  += psi_l_coef_bi_ortho(j,1) * hmono
+    accu_l_htwoe  += psi_l_coef_bi_ortho(j,1) * htwoe
+    accu_l_hthree += psi_l_coef_bi_ortho(j,1) * hthree
+    accu_l_htot   += psi_l_coef_bi_ortho(j,1) * htot
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+    accu_r_hmono  += psi_r_coef_bi_ortho(j,1) * hmono
+    accu_r_htwoe  += psi_r_coef_bi_ortho(j,1) * htwoe
+    accu_r_hthree += psi_r_coef_bi_ortho(j,1) * hthree
+    accu_r_htot   += psi_r_coef_bi_ortho(j,1) * htot
+   enddo
+   test_l_htot   += dabs(i_H_chi_array(0,1)-accu_l_htot)
+   test_l_hmono  += dabs(i_H_chi_array(1,1)-accu_l_hmono)
+   test_l_htwoe  += dabs(i_H_chi_array(2,1)-accu_l_htwoe)
+   test_l_hthree += dabs(i_H_chi_array(3,1)-accu_l_hthree)
+
+   test_r_htot   += dabs(i_H_phi_array(0,1)-accu_r_htot)
+   test_r_hmono  += dabs(i_H_phi_array(1,1)-accu_r_hmono)
+   test_r_htwoe  += dabs(i_H_phi_array(2,1)-accu_r_htwoe)
+   test_r_hthree += dabs(i_H_phi_array(3,1)-accu_r_hthree)
+
+ enddo
+ 
+ test_l_htot   *= 1.D0/dble(N_det)
+ test_l_hmono  *= 1.D0/dble(N_det)
+ test_l_htwoe  *= 1.D0/dble(N_det)
+ test_l_hthree *= 1.D0/dble(N_det)
+
+ test_r_htot   *= 1.D0/dble(N_det)
+ test_r_hmono  *= 1.D0/dble(N_det)
+ test_r_htwoe  *= 1.D0/dble(N_det)
+ test_r_hthree *= 1.D0/dble(N_det)
+ 
+ print*,'**************************'
+ print*,'test_l_htot    = ',test_l_htot
+ print*,'test_l_hmono   = ',test_l_hmono
+ print*,'test_l_htwoe   = ',test_l_htwoe
+ print*,'test_l_hthree  = ',test_l_hthree
+ print*,'**************************'
+ print*,'test_r_htot    = ',test_r_htot
+ print*,'test_r_hmono   = ',test_r_hmono
+ print*,'test_r_htwoe   = ',test_r_htwoe
+ print*,'test_r_hthree  = ',test_r_hthree
+
+end
+
+subroutine routine_grad_num
+ implicit none
+ integer :: indx,ihole,ipart
+ integer :: p,q
+ double precision :: accu_l, accu_r
+ double precision :: contrib_l, contrib_r
+
+ accu_l = 0.d0
+ accu_r = 0.d0
+  do indx=1,nMonoEx
+   q = excit(1,indx)
+   p = excit(2,indx)
+   contrib_l  =  dabs(dabs(gradvec_detail_left_old(0,indx)) - 2.D0 * dabs( Fock_matrix_tc_mo_tot(q,p)))
+   contrib_r  =  dabs(dabs(gradvec_detail_right_old(0,indx)) -2.D0 * dabs( Fock_matrix_tc_mo_tot(p,q)))
+   if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+    print*,indx,q,p
+    print*,gradvec_detail_left_old(0,indx),gradvec_detail_right_old(0,indx) 
+    print*,2.D0* Fock_matrix_tc_mo_tot(q,p), 2.d0* Fock_matrix_tc_mo_tot(p,q)
+   endif
+   accu_l += contrib_l
+   accu_r += contrib_r
+  enddo
+ print*,'accu_l,accu_r'
+ print*,accu_l,accu_r
+ 
+! do i = 1, nMonoEx
+!  print*,i,gradvec_old(i)
+! enddo
+
+end
+
+subroutine routine_grad_num_dm_one_body
+ implicit none
+ integer :: indx,ii,i,a,aa,tt,t,ibody
+ double precision :: accu_l, accu_r,ref_r, new_r, ref_l, new_l
+ double precision :: contrib_l, contrib_r
+ double precision :: res_l(0:3),res_r(0:3)
+
+ ibody = 2 ! check only the two-body term
+ provide gradvec_detail_left_old gradvec_tc_l 
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing inactive-->virtual'
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do ii = 1, n_core_inact_orb
+   do aa = 1, n_virt_orb
+    indx = mat_idx_c_v(ii,aa) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    i = list_core_inact(ii)
+    a = list_virt(aa)
+!    if(i==1.and.a==9)then
+!     print*,i,a,ref_r, new_r
+!     stop
+!    endif
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     print*,indx,i,a,ii,aa
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+     print*,gradvec_detail_left_old(0,indx),gradvec_tc_l(0,indx)
+     print*,gradvec_detail_right_old(0,indx),gradvec_tc_r(0,indx)
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+  print*,'**************************'
+  print*,'**************************'
+ endif
+ 
+ ibody = 2 ! check only the two-body term
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing inactive-->active'
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do ii = 1, n_core_inact_orb
+   do tt = 1, n_act_orb
+    indx = mat_idx_c_a(ii,tt) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     i = list_core_inact(ii)
+     t = list_act(tt)
+     print*,indx,i,t
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+ endif
+
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing active-->virtual '
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do tt = 1, n_act_orb
+   do aa = 1, n_virt_orb
+    indx = mat_idx_a_v(tt,aa) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     a = list_virt(aa)
+     t = list_act(tt)
+     print*,indx,t,a
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+!     print*,gradvec_detail_right_old(0,indx),gradvec_tc_r(0,indx)
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+ endif
+ 
+
+end

src/cisd/lccsd_prov.irp.f

@@ -3,7 +3,7 @@
  implicit none
  double precision, allocatable :: Dress_jj(:), H_jj(:), u_in(:,:)
  double precision :: ebefore, eafter, ecorr, thresh
- integer :: i,it
+ integer :: i,it,degree
  logical :: converged
  external H_u_0_nstates_openmp
  allocate(Dress_jj(N_det),H_jj(N_det),u_in(N_det,N_states_diag))
@@ -31,7 +31,11 @@
   print*,'ecorr = ',ecorr
   Dress_jj(1) = 0.d0
   do i = 2, N_det
-    if(ecorr + H_jj(i) .gt. H_jj(1))then
+    if(ecorr + H_jj(i) .lt. H_jj(1))then
+     print*,'Warning, some dets are not dressed: ' 
+     call get_excitation_degree(ref_bitmask,psi_det(1,1,i),degree,N_int)
+     print*,'degree, Delta E, coef', degree, H_jj(i)-H_jj(1), u_in(i,1)/u_in(1,1)
+    else
      Dress_jj(i) = ecorr
     endif
   enddo