Good Dyall Hamiltonian

config/ifort.cfg

@@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
 # 0 : Deactivate
 # 
 [OPTION]
-MODE    : DEBUG      ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
 CACHE   : 1          ; Enable cache_compile.py
 OPENMP  : 1          ; Append OpenMP flags
 

plugins/Full_CI/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_full
+Perturbation Selectors_full Generators_full MRPT_Utils

plugins/Full_CI/micro_pt2.irp.f

@@ -1,42 +0,0 @@
-program micro_pt2
-  implicit none
-  BEGIN_DOC
-! Helper program to compute the PT2 in distributed mode.
-  END_DOC
-
-  read_wf = .False.
-  SOFT_TOUCH read_wf
-  call provide_everything
-  call switch_qp_run_to_master
-  call run_wf
-
-end
-
-subroutine provide_everything
-  PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
-end
-
-subroutine run_wf
-  use f77_zmq
-  implicit none
-
-  integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
-  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
-
-  print *,  'Getting wave function'
-  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
-
-  call zmq_get_psi(zmq_to_qp_run_socket, 1)
-  call write_double(6,ci_energy,'Energy')
-  zmq_state = 'h_apply_fci_pt2'
-
-  call provide_everything
-  integer :: rc, i
-
-  !$OMP PARALLEL PRIVATE(i)
-  i = omp_get_thread_num()
-  call H_apply_FCI_PT2_slave_tcp(i)
-  !$OMP END PARALLEL
-
-
-end

plugins/MRPT_Utils/MRPT_Utils.main.irp.f

@@ -16,8 +16,8 @@ subroutine routine_3
  print *,  'N_det    = ', N_det
  print *,  'N_states = ', N_states
  print *,  'PT2      = ', second_order_pt_new(1)
- print *,  'E        = ', CI_energy
+ print *,  'E        = ', CI_energy(1)
- print *,  'E+PT2    = ', CI_energy+second_order_pt_new(1)
+ print *,  'E+PT2    = ', CI_energy(1)+second_order_pt_new(1)
  print *,'****** DIAGONALIZATION OF DRESSED MATRIX ******'
  print *,  'E dressed= ', CI_dressed_pt2_new_energy(1)
 

plugins/MRPT_Utils/energies_cas.irp.f

@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
  integer :: i 
  double precision :: energies(N_states_diag)
  do i = 1, N_states
-  call u0_H_dyall_u0(energies,psi_det,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
+  call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
   energy_cas_dyall(i) = energies(i)
   print*,  'energy_cas_dyall(i)',  energy_cas_dyall(i)
  enddo
@@ -33,8 +33,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
       psi_in_out_coef(i,j) = psi_coef(i,j)
     enddo
     do j = 1, N_int
-     psi_in_out(j,1,i) =  psi_det(j,1,i) 
+     psi_in_out(j,1,i) =  psi_active(j,1,i) 
-     psi_in_out(j,2,i) =  psi_det(j,2,i) 
+     psi_in_out(j,2,i) =  psi_active(j,2,i) 
     enddo
    enddo
     do  state_target = 1,N_states
@@ -71,8 +71,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
       psi_in_out_coef(i,j) = psi_coef(i,j)
     enddo
     do j = 1, N_int
-     psi_in_out(j,1,i) =  psi_det(j,1,i) 
+     psi_in_out(j,1,i) =  psi_active(j,1,i) 
-     psi_in_out(j,2,i) =  psi_det(j,2,i) 
+     psi_in_out(j,2,i) =  psi_active(j,2,i) 
     enddo
    enddo
    do state_target = 1, N_states
@@ -115,8 +115,8 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
         psi_in_out_coef(i,j) = psi_coef(i,j)
       enddo
       do j = 1, N_int
-       psi_in_out(j,1,i) =  psi_det(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
-       psi_in_out(j,2,i) =  psi_det(j,2,i) 
+       psi_in_out(j,2,i) =  psi_active(j,2,i) 
       enddo
      enddo
      do state_target = 1 , N_states
@@ -164,8 +164,8 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
         psi_in_out_coef(i,j) = psi_coef(i,j)
       enddo
       do j = 1, N_int
-       psi_in_out(j,1,i) =  psi_det(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
-       psi_in_out(j,2,i) =  psi_det(j,2,i) 
+       psi_in_out(j,2,i) =  psi_active(j,2,i) 
       enddo
      enddo
      call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
@@ -210,8 +210,8 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
         psi_in_out_coef(i,j) = psi_coef(i,j)
       enddo
       do j = 1, N_int
-       psi_in_out(j,1,i) =  psi_det(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
-       psi_in_out(j,2,i) =  psi_det(j,2,i) 
+       psi_in_out(j,2,i) =  psi_active(j,2,i) 
       enddo
      enddo
      do state_target = 1, N_states
@@ -265,8 +265,8 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
           psi_in_out_coef(i,j) = psi_coef(i,j)
         enddo
         do j = 1, N_int
-         psi_in_out(j,1,i) =  psi_det(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
-         psi_in_out(j,2,i) =  psi_det(j,2,i) 
+         psi_in_out(j,2,i) =  psi_active(j,2,i) 
         enddo
        enddo
 
@@ -325,8 +325,8 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
           psi_in_out_coef(i,j) = psi_coef(i,j)
         enddo
         do j = 1, N_int
-         psi_in_out(j,1,i) =  psi_det(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
-         psi_in_out(j,2,i) =  psi_det(j,2,i) 
+         psi_in_out(j,2,i) =  psi_active(j,2,i) 
         enddo
        enddo
        do state_target = 1, N_states
@@ -384,8 +384,8 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
           psi_in_out_coef(i,j) = psi_coef(i,j)
         enddo
         do j = 1, N_int
-         psi_in_out(j,1,i) =  psi_det(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
-         psi_in_out(j,2,i) =  psi_det(j,2,i) 
+         psi_in_out(j,2,i) =  psi_active(j,2,i) 
         enddo
        enddo
        do state_target = 1, N_states
@@ -443,8 +443,8 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
           psi_in_out_coef(i,j) = psi_coef(i,j)
         enddo
         do j = 1, N_int
-         psi_in_out(j,1,i) =  psi_det(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
-         psi_in_out(j,2,i) =  psi_det(j,2,i) 
+         psi_in_out(j,2,i) =  psi_active(j,2,i) 
         enddo
        enddo
        do state_target = 1, N_states

plugins/MRPT_Utils/excitations_cas.irp.f

@@ -130,6 +130,24 @@ double precision function diag_H_mat_elem_no_elec_check(det_in,Nint)
    enddo
   enddo 
 
+  ! alpha - core-act
+  do i = 1, elec_num_tab_local(1)
+   iorb =  occ(i,1)
+   do j = 1, n_core_inact_orb
+    jorb = list_core_inact(j)
+    diag_H_mat_elem_no_elec_check +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) - mo_bielec_integral_jj_exchange(jorb,iorb)
+   enddo
+  enddo 
+
+  ! beta - core-act
+  do i = 1, elec_num_tab_local(2)
+   iorb =  occ(i,2)
+   do j = 1, n_core_inact_orb
+    jorb = list_core_inact(j)
+    diag_H_mat_elem_no_elec_check +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) - mo_bielec_integral_jj_exchange(jorb,iorb)
+   enddo
+  enddo 
+  
 end
 
 subroutine a_operator_no_check(iorb,ispin,key,hjj,Nint,na,nb)

plugins/Perturbation/pt2_new.irp.f

@@ -42,17 +42,7 @@ subroutine i_H_psi_pert_new_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,
       i_H_psi_array(1) = i_H_psi_array(1) + coef(i_in_coef,1)*hij
       call get_delta_e_dyall(keys(1,1,i_in_key),key,delta_e_final)
 
-      if(delta_e_final == 0.d0)then
-       call get_delta_e_dyall_verbose(keys(1,1,i_in_key),key,delta_e_final)
-       call debug_det(keys(1,1,i_in_key),N_int)
-       call debug_det(key,N_int)
-       stop
-      endif
       coef_pert +=  coef(i_in_coef,1)*hij / delta_e_final
-!     print*, 'delta_e_final = ',delta_e_final
-
-!     call i_H_j(key,key,Nint,hjj)
-!     coef_pert +=  coef(i_in_coef,1)*hij / (CI_electronic_energy(1) - hjj)
     enddo
     if     (coef_pert * i_H_psi_array(1) > 0.d0)then
       print*, coef_pert * i_H_psi_array(1)

plugins/loc_cele/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-MO_Basis
+MO_Basis  Integrals_Bielec Bitmask

plugins/loc_cele/loc_cele.irp.f

@@ -92,7 +92,7 @@
 
 
 
-      nrot(1) = 6   ! number of orbitals to be localized
+      nrot(1) = 2   ! number of orbitals to be localized
 
 
       integer :: index_rot(1000,1)
@@ -101,9 +101,10 @@
        cmoref = 0.d0
        irot = 0
 
-      do i=1,nrot(1)
+       irot(1,1) = 11
-       irot(i,1) = 19+i 
+       irot(2,1) = 12
-      enddo
+       cmoref(15,1,1)   = 1.d0   ! 
+       cmoref(14,2,1)   = 1.d0   ! 
 
 ! ESATRIENE with 3 bonding and anti bonding orbitals 
 ! First bonding orbital for esa

src/Determinants/diagonalize_CI.irp.f

@@ -91,6 +91,8 @@ END_PROVIDER
       do j=1,N_det
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
         s2_eigvalues(j) = s2
+        print*, 's2 in lapack',s2
+        print*, eigenvalues(j)
         ! Select at least n_states states with S^2 values closed to "expected_s2"
         if(dabs(s2-expected_s2).le.0.3d0)then
          i_state +=1