minor modif

src/fci/fci.irp.f

@@ -42,7 +42,7 @@ program fci
     write(json_unit,json_array_open_fmt) 'fci'
 
     double precision, allocatable :: Ev(:),PT2(:)
-    allocate(Ev(N_states), PT2(N_state))
+    allocate(Ev(N_states), PT2(N_states))
     if (do_pt2) then
       call run_stochastic_cipsi(Ev,PT2)
     else

src/non_hermit_dav/biorthog.irp.f

@@ -270,7 +270,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   integer,          intent(out) :: n_real_eigv
   double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
 
-  integer                       :: i, j
+  integer                       :: i, j,k 
   integer                       :: n_good
   double precision              :: thr, thr_cut, thr_diag, thr_norm
   double precision              :: accu_d, accu_nd
@@ -278,6 +278,8 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   integer,          allocatable :: list_good(:), iorder(:)
   double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
   double precision, allocatable :: S(:,:)
+  double precision, allocatable  :: phi_1_tilde(:),phi_2_tilde(:),chi_1_tilde(:),chi_2_tilde(:)
+  allocate(phi_1_tilde(n),phi_2_tilde(n),chi_1_tilde(n),chi_2_tilde(n))
 
 
   ! -------------------------------------------------------------------------------------
@@ -301,11 +303,78 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   call lapack_diag_non_sym(n, A, WR, WI, VL, VR)
   !call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)
 
-  !print *, ' '
+
-  !print *, ' eigenvalues'
+
-  !do i = 1, n
+  print *, ' '
-  !  write(*, '(1000(F16.10,X))') WR(i), WI(i)
+  print *, ' eigenvalues'
-  !enddo
+  i = 1
+  do while(i .le. n)
+    write(*, '(I3,X,1000(F16.10,X))')i, WR(i), WI(i)
+   if(.false.)then
+    if(WI(i).ne.0.d0)then
+     print*,'*****************'
+     print*,'WARNING ! IMAGINARY EIGENVALUES !!!'
+     write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
+     ! phi = VR(:,i), psi = VR(:,i+1), |Phi_i> = phi + j psi , |Phi_i+1> = phi - j psi 
+     ! chi = VL(:,i), xhi = VL(:,i+1), |Chi_i> = chi + j xhi , |Chi_i+1> = chi - j xhi 
+     !                                   
+     accu_chi_phi = 0.d0
+     accu_xhi_psi = 0.d0
+     accu_chi_psi = 0.d0
+     accu_xhi_phi = 0.d0
+     double precision :: accu_chi_phi, accu_xhi_psi, accu_chi_psi, accu_xhi_phi
+     double precision :: mat_ovlp(2,2),eigval_tmp(2),eigvec(2,2),mat_ovlp_orig(2,2)
+     do j = 1, n 
+      accu_chi_phi += VL(j,i)   * VR(j,i)
+      accu_xhi_psi += VL(j,i+1) * VR(j,i+1)
+      accu_chi_psi += VL(j,i)   * VR(j,i+1)
+      accu_xhi_phi += VL(j,i+1) * VR(j,i)
+     enddo
+     mat_ovlp_orig(1,1) = accu_chi_phi
+     mat_ovlp_orig(2,1) = accu_xhi_phi
+     mat_ovlp_orig(1,2) = accu_chi_psi
+     mat_ovlp_orig(2,2) = accu_xhi_psi
+     print*,'old overlap matrix '
+     write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,1)
+     write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,2)
+
+
+     mat_ovlp(1,1) = accu_xhi_phi
+     mat_ovlp(2,1) = accu_chi_phi
+     mat_ovlp(1,2) = accu_xhi_psi
+     mat_ovlp(2,2) = accu_chi_psi
+    !print*,'accu_chi_phi = ',accu_chi_phi
+    !print*,'accu_xhi_psi = ',accu_xhi_psi
+    !print*,'accu_chi_psi = ',accu_chi_psi
+    !print*,'accu_xhi_phi = ',accu_xhi_phi
+     print*,'new overlap matrix '
+     write(*,'(100(F16.10,X))')mat_ovlp(1:2,1)
+     write(*,'(100(F16.10,X))')mat_ovlp(1:2,2)
+     call lapack_diag(eigval_tmp,eigvec,mat_ovlp,2,2)
+     print*,'eigval_tmp(1) = ',eigval_tmp(1)
+     print*,'eigvec(1) = ',eigvec(1:2,1)
+     print*,'eigval_tmp(2) = ',eigval_tmp(2)
+     print*,'eigvec(2) = ',eigvec(1:2,2)
+     print*,'*****************'
+     phi_1_tilde = 0.d0
+     phi_2_tilde = 0.d0
+     chi_1_tilde = 0.d0
+     chi_2_tilde = 0.d0
+     do j = 1, n
+      phi_1_tilde(j) += VR(j,i) * eigvec(1,1) + VR(j,i+1) * eigvec(2,1)
+      phi_2_tilde(j) += VR(j,i) * eigvec(1,2) + VR(j,i+1) * eigvec(2,2)
+      chi_1_tilde(j) += VL(j,i+1) * eigvec(1,1) + VL(j,i) * eigvec(2,1)
+      chi_2_tilde(j) += VL(j,i+1) * eigvec(1,2) + VL(j,i) * eigvec(2,2)
+     enddo
+     VR(1:n,i)   = phi_1_tilde(1:n)
+     VR(1:n,i+1) = phi_2_tilde(1:n)
+!     Vl(1:n,i)   = -chi_1_tilde(1:n)
+!     Vl(1:n,i+1) = chi_2_tilde(1:n)
+     i+=1 
+    endif
+   endif
+    i+=1
+  enddo
   !print *, ' right eigenvect bef' 
   !do i = 1, n
   !  write(*, '(1000(F16.10,X))') VR(:,i)
@@ -429,6 +498,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
 !   call impose_orthog_degen_eigvec(n, eigval, reigvec)
 !   call impose_orthog_degen_eigvec(n, eigval, leigvec)
 
+    call reorder_degen_eigvec(n, eigval, leigvec, reigvec)
     call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
 
 

src/non_hermit_dav/lapack_diag_non_hermit.irp.f

@@ -1944,6 +1944,96 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S)
 end subroutine check_orthog
 
 ! ---
+subroutine reorder_degen_eigvec(n, e0, L0, R0)
+
+  implicit none
+
+  integer,          intent(in)    :: n
+  double precision, intent(in)    :: e0(n)
+  double precision, intent(inout) :: L0(n,n), R0(n,n)
+
+  logical                         :: complex_root
+  integer                         :: i, j, k, m
+  double precision                :: ei, ej, de, de_thr
+  double precision                :: accu_d, accu_nd
+  integer,          allocatable   :: deg_num(:)
+  double precision, allocatable   :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
+
+  ! ---
+
+  allocate( deg_num(n) )
+  do i = 1, n
+    deg_num(i) = 1
+  enddo
+
+  de_thr = thr_degen_tc
+
+  do i = 1, n-1
+    ei = e0(i)
+
+    ! already considered in degen vectors
+    if(deg_num(i).eq.0) cycle
+
+    do j = i+1, n
+      ej = e0(j)
+      de = dabs(ei - ej)
+
+      if(de .lt. de_thr) then
+        deg_num(i) = deg_num(i) + 1 
+        deg_num(j) = 0
+      endif
+
+    enddo
+  enddo
+  
+  do i = 1, n
+    if(deg_num(i) .gt. 1) then
+      print *, ' degen on', i, deg_num(i), e0(i)
+    endif
+  enddo
+
+  ! ---
+
+  do i = 1, n
+    m = deg_num(i)
+
+    if(m .gt. 1) then
+  
+      allocate(L(n,m))
+      allocate(R(n,m),S(m,m))
+
+      do j = 1, m
+        L(1:n,j) = L0(1:n,i+j-1)
+        R(1:n,j) = R0(1:n,i+j-1)
+      enddo
+
+      call dgemm( 'T', 'N', m, m, n, 1.d0      &
+                , L, size(L, 1), R, size(R, 1) &
+                , 0.d0, S, size(S, 1) )
+      print*,'Overlap matrix '
+      accu_nd = 0.D0
+      do j = 1, m
+       write(*,'(100(F16.10,X))')S(1:m,j)
+       do k = 1, m
+        if(j==k)cycle
+        accu_nd += dabs(S(j,k))
+       enddo
+      enddo
+      print*,'accu_nd = ',accu_nd
+!      if(accu_nd .gt.1.d-10)then
+!        stop
+!      endif
+      do j = 1, m
+        L0(1:n,i+j-1) = L(1:n,j)
+        R0(1:n,i+j-1) = R(1:n,j)
+      enddo
+
+      deallocate(L, R,S)
+
+    endif
+  enddo
+
+end subroutine reorder_degen_eigvec 
 
 subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
@@ -2030,7 +2120,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
       !call bi_ortho_s_inv_half(m, L, R, S_inv_half)
       !deallocate(S, S_inv_half)
 
-!      call impose_biorthog_svd(n, m, L, R)
+      call impose_biorthog_svd(n, m, L, R)
 !     call impose_biorthog_inverse(n, m, L, R)
 
       !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
@@ -2046,7 +2136,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
     endif
   enddo
-  call impose_biorthog_inverse(n, n, L0, R0)
+!  call impose_biorthog_inverse(n, n, L0, R0)
 
 end subroutine impose_biorthog_degen_eigvec 
 

src/tc_bi_ortho/tc_natorb.irp.f

@@ -32,6 +32,13 @@
   thr_d   = 1.d-6
   thr_nd  = 1.d-6
   thr_deg = 1.d-3
+  do i = 1, mo_num
+   do j = 1, mo_num
+    if(dabs(dm_tmp(j,i)).lt.thr_d)then
+     dm_tmp(j,i) = 0.d0
+    endif
+   enddo
+  enddo
 !  if(n_core_orb.ne.0)then
 !   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
 !                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)