Merge pull request #6 from QuantumPackage/dev

Changes Manus

bin/qp_reset

@@ -105,6 +105,7 @@ if [[ $mos -eq 1 ]] ; then
       echo "Warning: You will need to re-define the MO classes"
   fi
   rm --recursive --force -- ${ezfio}/mo_basis
+  rm --recursive --force -- ${ezfio}/bi_ortho_mos
   rm --recursive --force -- ${ezfio}/work/mo_ints_*
 fi
 

external/qp2-dependencies

@@ -1 +1 @@
-Subproject commit 242151e03d1d6bf042387226431d82d35845686a
+Subproject commit 9e5b27ce5a174901765cec9db9e7b2aa6170a5de

src/ao_many_one_e_ints/list_grid.irp.f

@@ -8,8 +8,8 @@
  thr = 1.d-11
  print*,' expo_good_j_mu_1gauss = ',expo_good_j_mu_1gauss
  !$OMP PARALLEL DEFAULT (NONE)                                      &
- !$OMP PRIVATE (ipoint, i, j, r, overlap, thr,fact_gauss, alpha, center,dist,sigma,center_ij) &
+ !$OMP PRIVATE (ipoint, i, j, r, overlap, fact_gauss, alpha, center,dist,sigma,center_ij) &
- !$OMP SHARED  (n_points_final_grid, ao_num, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,&
+ !$OMP SHARED  (n_points_final_grid, ao_num, thr, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,&
  !$OMP          final_grid_points,ao_prod_center,ao_prod_sigma,ao_nucl)
  !$OMP DO
  do i = 1, ao_num

src/cipsi_tc_bi_ortho/selection.irp.f

@@ -771,12 +771,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         ! <det|H(j)|psi_0> and transpose
         ! -------------------------------------------
 
-        istate = 1
 !        call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii)
         double precision               :: hmono, htwoe, hthree
         call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
+        do istate = 1,N_states
           delta_E = E0(istate) - Hii + E_shift
-
           double precision               :: alpha_h_psi_tmp, psi_h_alpha_tmp, error
           if(debug_tc_pt2 == 1)then !! Using the old version
             psi_h_alpha = 0.d0
@@ -797,85 +796,37 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
               call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
               psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
               alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
-!            psi_h_alpha += i_h_alpha * 1.d0 ! left function 
-!            alpha_h_psi += alpha_h_i * 1.d0 ! right function
             enddo
             if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then
               error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi)
               if(error.gt.1.d-2)then
                 print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E
               endif
-!           if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-08 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-08)then
-!               call debug_det(det,N_int)
-!               print*,'psi_h_alpha,alpha_h_psi'
-!               print*,psi_h_alpha,alpha_h_psi
-!               print*,psi_h_alpha_tmp,alpha_h_psi_tmp
-!               print*,dabs(psi_h_alpha - psi_h_alpha_tmp),dabs(alpha_h_psi - alpha_h_psi_tmp)
-!               alpha_h_psi = 0.d0
-!               psi_h_alpha = 0.d0
-!               do iii = 1, N_det
-!                 
-!                 call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int)
-!                 call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha)
-!                 call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i)
-!                 alpha_h_psi += alpha_h_i
-!                 psi_h_alpha += i_h_alpha
-!                 if(dabs(i_h_alpha).gt.1.d-10.or.dabs(alpha_h_i).gt.1.d-10)then
-!                  call debug_det(psi_det(1,1,iii),N_int)
-!                  print*,iii,degree,i_h_alpha,alpha_h_i
-!                  print*,psi_h_alpha,alpha_h_psi
-!                  print*,leigvec_tc_bi_orth(iii,1),reigvec_tc_bi_orth(iii,1)
-!                 endif
-!               enddo
-!               stop
-!            endif
             endif
           else
            psi_h_alpha = mat_l(istate, p1, p2)
            alpha_h_psi = mat_r(istate, p1, p2)
           endif
-
-        !if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then
-        !  print*, 'E0,Hii,E_shift'
-        !  print*, E0(istate), Hii, E_shift
-        !  print*, psi_h_alpha, alpha_h_psi, delta_E
-        !  print*, psi_h_alpha * alpha_h_psi / delta_E
-        !  !if(Hii .lt. E0(istate)) then
-        !  !  call debug_det(det, N_int)
-        !  !  print*, ' |E0| < |Hii| !!!'
-        !  !  print*, ' E0  = ', E0(istate)
-        !  !  print*, ' Hii = ', Hii
-        !  !endif
-        !endif
-
           coef(istate)   = alpha_h_psi / delta_E 
           e_pert(istate) = coef(istate) * psi_h_alpha
-        if(selection_tc     ==  1 )then
+!         if(selection_tc     ==  1 )then
-         if(e_pert(istate).lt.0.d0)then
+!          if(e_pert(istate).lt.0.d0)then
-          e_pert(istate) = 0.d0
+!           e_pert(istate) = 0.d0
-         endif
-        else if(selection_tc == -1)then
-         if(e_pert(istate).gt.0.d0)then
-          e_pert(istate) = 0.d0
-         endif
-        endif
-         
-
-        !if(e_pert(istate) .gt. 1.d-15) then
-        !  print*, 'E0,Hii,E_shift'
-        !  print*, E0(istate), Hii, E_shift
-        !  print*, psi_h_alpha, alpha_h_psi, delta_E
-        !  print*, psi_h_alpha*alpha_h_psi/delta_E
 !          endif
-
+!         else if(selection_tc == -1)then
-!      elseif(cipsi_tc == "h_tc_2x2") then
+!          if(e_pert(istate).gt.0.d0)then
+!           e_pert(istate) = 0.d0
+!          endif
+!         endif
+        enddo
 
 
       do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1
 
       do istate = 1, N_states
 
-        alpha_h_psi = mat(istate, p1, p2)
+        alpha_h_psi = mat_r(istate, p1, p2)
+        psi_h_alpha = mat_l(istate, p1, p2)
 
         pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate)
         pt2_data % variance(istate)  = pt2_data % variance(istate) + dabs(e_pert(istate))
@@ -885,10 +836,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           case(5)
             ! Variance selection
             if (h0_type == 'CFG') then
-              w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) & 
+              w = min(w, - psi_h_alpha * alpha_h_psi * s_weight(istate,istate)) &
                 / c0_weight(istate)
             else
-              w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate))
+              w = min(w, - psi_h_alpha * alpha_h_psi * s_weight(istate,istate))
             endif
           case(6)
             if (h0_type == 'CFG') then

src/csf/sigma_vector.irp.f

@@ -236,7 +236,7 @@ end subroutine get_phase_qp_to_cfg
   ! initialization
   psi_coef_config = 0.d0
   DetToCSFTransformationMatrix(0,:,:) = 1.d0
-  do i = 2-iand(MS,1_8), NSOMOMax,2
+  do i = 2-int(iand(MS,1_8),4), NSOMOMax,2
     Isomo = IBSET(0_8, i) - 1_8
     ! rows = Ncsfs
     ! cols = Ndets

src/fci_tc_bi/copy_wf.irp.f

@@ -100,7 +100,7 @@ subroutine copy_H_apply_buffer_to_wf_tc
 
   logical :: found_duplicates
   call remove_duplicates_in_psi_det_tc(found_duplicates)
-  call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,size(psi_l_coef_bi_ortho,1),N_det,N_states)
+  call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_det,size(psi_l_coef_bi_ortho,1),N_states)
   SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho
 
 end

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -100,33 +100,33 @@
     print*,'Computing the left-eigenvector '
     vec_tmp = 0.d0
     do istate = 1, N_states
-     vec_tmp(:,istate) = psi_l_coef_bi_ortho(:,istate)
+     vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
     enddo
     do istate = N_states+1, n_states_diag
      vec_tmp(istate,istate) = 1.d0
     enddo
     call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
     do istate = 1, N_states
-     leigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
+     leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo
 
     print*,'Computing the right-eigenvector '
  !!!! Preparing the right-eigenvector
     vec_tmp = 0.d0
     do istate = 1, N_states
-     vec_tmp(:,istate) = psi_r_coef_bi_ortho(:,istate)
+     vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
     enddo
     do istate = N_states+1, n_states_diag
      vec_tmp(istate,istate) = 1.d0
     enddo
     call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
     do istate = 1, N_states
-     reigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
+     reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo
 
     deallocate(H_jj)
    endif
-  call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,N_det,N_det,N_states)
+  call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states)
    print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
    norm_ground_left_right_bi_orth = 0.d0
    do j = 1, N_det

src/tc_bi_ortho/tc_natorb.irp.f

@@ -14,7 +14,7 @@
   END_DOC
 
   implicit none
-  integer                       :: i, j, k, n_real
+  integer                       :: i, j, k
   double precision              :: thr_d, thr_nd, thr_deg, accu
   double precision              :: accu_d, accu_nd
   double precision, allocatable :: dm_tmp(:,:), fock_diag(:)
@@ -36,17 +36,17 @@
                               , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
 !   call non_hrmt_bieig( mo_num, dm_tmp&
 !                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
-!                      , n_real, natorb_tc_eigval )
+!                      , mo_num, natorb_tc_eigval )
 
   accu = 0.d0
-  do i = 1, n_real
+  do i = 1, mo_num
     print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
     accu += -natorb_tc_eigval(i)
   enddo
   print *, ' accu = ', accu
 
   dm_tmp = 0.d0
-  do i = 1, n_real
+  do i = 1, mo_num
     accu = 0.d0
     do k = 1, mo_num
       accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
@@ -54,7 +54,7 @@
     accu = 1.d0/dsqrt(dabs(accu))
     natorb_tc_reigvec_mo(:,i) *= accu
     natorb_tc_leigvec_mo(:,i) *= accu
-    do j = 1, n_real
+    do j = 1, mo_num
       do k = 1, mo_num
         dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
       enddo

src/tc_keywords/EZFIO.cfg

@@ -94,12 +94,6 @@ doc: Maximum number of SCF iterations
 interface: ezfio,provider,ocaml
 default: 100
 
-[selection_tc]
-type: integer
-doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative
-interface: ezfio,provider,ocaml
-default: 0
-
 [j1b_pen]
 type: double precision
 doc: exponents of the 1-body Jastrow