pouet

plugins/MRPT/MRPT_Utils.main.irp.f

@@ -19,8 +19,8 @@ subroutine routine_3
  do i = 1, N_States
   print*,'State',i
   write(*,'(A12,X,I3,A3,XX,F16.10)')    '  PT2       ', i,' = ', second_order_pt_new(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E         ', i,' = ', CI_energy(i)
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E         ', i,' = ', psi_ref_average_value(i)
-  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E+PT2     ', i,' = ', CI_energy(i)+second_order_pt_new(i)
+  write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E+PT2     ', i,' = ', psi_ref_average_value(i)+second_order_pt_new(i)
   write(*,'(A12,X,I3,A3,XX,F16.09)')    '  E dressed ', i,' = ', CI_dressed_pt2_new_energy(i)
   write(*,'(A12,X,I3,A3,XX,F16.09)')    '  S^2       ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
   print*,'coef before and after'
@@ -28,6 +28,11 @@ subroutine routine_3
    print*,psi_ref_coef(j,i),CI_dressed_pt2_new_eigenvectors(j,i)
   enddo
  enddo 
+ if(save_heff_eigenvectors)then
+  call save_wavefunction_general(N_det_ref,N_states_diag_heff,psi_ref,N_det_ref,CI_dressed_pt2_new_eigenvectors)
+ endif
+! print*, 'neutral = ',psi_ref_coef(1,1),CI_dressed_pt2_new_eigenvectors(1,1)
+! print*, 'ionic   = ',psi_ref_coef(3,1),CI_dressed_pt2_new_eigenvectors(3,1)
 
 end
 

plugins/MRPT/print_1h2p.irp.f

@@ -9,11 +9,19 @@ subroutine routine_2
  implicit none
  integer :: i,j,degree
  double precision :: hij
-!provide   one_creat_virt
+ do i =1, n_core_inact_orb 
- do i =1, n_act_orb
+  write(*,'(I3,x,100(F16.10,X))')list_core_inact(i),fock_core_inactive_total_spin_trace(list_core_inact(i),1)
-  write(*,'(I3,x,100(F16.10,X))')i,one_creat(i,:,1)
+ enddo
-! write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(1,4,1,2,1)
+ print*,''
-! 
+ do i =1, n_virt_orb
+  write(*,'(I3,x,100(F16.10,X))')list_virt(i),fock_virt_total_spin_trace(list_virt(i),1)
+ enddo
+ stop
+ do i = 1, n_virt_orb
+  do j = 1, n_inact_orb 
+  if(dabs(one_anhil_one_creat_inact_virt(j,i,1)) .lt. 1.d-10)cycle
+  write(*,'(I3,x,I3,X,100(F16.10,X))')list_virt(i),list_inact(j),one_anhil_one_creat_inact_virt(j,i,1)
+  enddo
  enddo
 
 

plugins/MRPT_Utils/EZFIO.cfg

@@ -5,6 +5,13 @@ interface: ezfio,provider,ocaml
 default: True
 
 
+[save_heff_eigenvectors]
+type: logical
+doc: If true, you save the eigenvectors of the effective hamiltonian
+interface: ezfio,provider,ocaml
+default: False
+
+
 [pure_state_specific_mrpt2]
 type: logical
 doc: If true, diagonalize the dressed matrix for each state and do a state following of the initial states
@@ -12,3 +19,9 @@ interface: ezfio,provider,ocaml
 default: True
 
 
+[N_states_diag_heff]
+type: States_number
+doc: Number of eigenvectors obtained with the effective hamiltonian
+interface: ezfio,provider,ocaml
+default: 1
+

plugins/MRPT_Utils/energies_cas.irp.f

@@ -617,6 +617,9 @@ END_PROVIDER
   
  thresh_norm = 1.d-20
 
+!do i = 1, N_det_ref
+! print*, psi_ref_coef(i,1)
+!enddo
 
 
  do vorb = 1,n_virt_orb
@@ -645,6 +648,10 @@ END_PROVIDER
         double precision ::  coef,contrib
         coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j)
         psi_in_out_coef(i,j) = coef * hij 
+!       if(vorb == 1.and. iorb == 1)then
+!       if(vorb == 1.and. iorb == 3)then
+!       print*, i,hij,coef
+!       endif
         norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) 
       enddo
      enddo

plugins/MRPT_Utils/ezfio_interface.irp.f

@@ -22,6 +22,44 @@ BEGIN_PROVIDER [ logical, do_third_order_1h1p  ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ logical, save_heff_eigenvectors  ]
+  implicit none
+  BEGIN_DOC
+! If true, you save the eigenvectors of the effective hamiltonian
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrpt_utils_save_heff_eigenvectors(has)
+  if (has) then
+    call ezfio_get_mrpt_utils_save_heff_eigenvectors(save_heff_eigenvectors)
+  else
+    print *, 'mrpt_utils/save_heff_eigenvectors not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_states_diag_heff  ]
+  implicit none
+  BEGIN_DOC
+! Number of eigenvectors obtained with the effective hamiltonian
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrpt_utils_n_states_diag_heff(has)
+  if (has) then
+    call ezfio_get_mrpt_utils_n_states_diag_heff(n_states_diag_heff)
+  else
+    print *, 'mrpt_utils/n_states_diag_heff not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
 BEGIN_PROVIDER [ logical, pure_state_specific_mrpt2  ]
   implicit none
   BEGIN_DOC

plugins/MRPT_Utils/mrpt_dress.irp.f

@@ -117,19 +117,15 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
        do i_state = 1, N_states
         delta_e(i_state) = 1.d+20
        enddo
-     !else if(degree_scalar== 1)then
       else 
        call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
-      !if(dabs(delta_e(2)) .le. dabs(0.01d0))then
-      !print*, delta_e(2)
-      !call debug_det(psi_ref(1,1,index_i),N_int)
-      !call debug_det(tq(1,1,i_alpha),N_int)
-      !endif
       
-      !else 
+     ! !!!!!!!!!!!!! SHIFTED BK 
-      !do i_state = 1, N_states
+     ! double precision :: hjj
-      ! delta_e(i_state) = 1.d+20
+     ! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
-      !enddo
+     ! delta_e(1) = CI_electronic_energy(1) - hjj
+     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       endif
       hij_array(index_i) = hialpha
       do i_state = 1,N_states

plugins/MRPT_Utils/mrpt_utils.irp.f

@@ -48,6 +48,8 @@
    do i = 1, N_det_ref
      write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state)
     do j = 1, N_det_ref
+!    print*, accu
+!    print*,delta_ij_tmp(j,i,i_state) , psi_ref_coef(i,i_state) , psi_ref_coef(j,i_state)
      accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
      delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
     enddo
@@ -65,11 +67,41 @@
      write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state)
     do j = 1, N_det_ref
      accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
-      delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+     delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+    enddo
+   enddo
+   double precision :: accu_diag,accu_non_diag
+   accu_diag = 0.d0 
+   accu_non_diag = 0.d0
+   do i = 1, N_det_ref
+    accu_diag += delta_ij_tmp(i,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(i,i_state)
+    do j = 1, N_det_ref
+     if(i == j)cycle
+     accu_non_diag += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
     enddo 
    enddo
    second_order_pt_new_1h1p(i_state) = accu(i_state) 
    enddo
+  !double precision :: neutral, ionic
+  !neutral = 0.d0
+  !do i = 1, 2 
+  ! do j = 1, N_det_ref
+  !  neutral += psi_ref_coef(j,1) * delta_ij_tmp(j,i,1) * psi_ref_coef(i,1)
+  ! enddo
+  !enddo
+  !do i = 3, 4 
+  ! do j = 1, N_det_ref
+  !  ionic += psi_ref_coef(j,1) * delta_ij_tmp(j,i,1) * psi_ref_coef(i,1)
+  ! enddo
+  !enddo
+  !neutral = delta_ij_tmp(1,1,1) * psi_ref_coef(1,1)**2 + delta_ij_tmp(2,2,1) * psi_ref_coef(2,1)**2 & 
+  !        + delta_ij_tmp(1,2,1) * psi_ref_coef(1,1)* psi_ref_coef(2,1) + delta_ij_tmp(2,1,1) * psi_ref_coef(1,1)* psi_ref_coef(2,1)
+  !ionic   = delta_ij_tmp(3,3,1) * psi_ref_coef(3,1)**2 + delta_ij_tmp(4,4,1) * psi_ref_coef(4,1)**2 & 
+  !        + delta_ij_tmp(3,4,1) * psi_ref_coef(3,1)* psi_ref_coef(4,1) + delta_ij_tmp(4,3,1) * psi_ref_coef(3,1)* psi_ref_coef(4,1)
+  !neutral = delta_ij_tmp(1,1,1) 
+  !ionic   = delta_ij_tmp(3,3,1) 
+  !print*, 'neutral = ',neutral
+  !print*, 'ionic   = ',ionic
    print*, '1h1p = ',accu
   
    ! 1h1p third order
@@ -168,6 +200,22 @@
    enddo
    print*, '2h2p = ',contrib_2h2p(:) 
 
+  !! 2h2p   old fashion
+  !delta_ij_tmp = 0.d0
+  !call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref)
+  !accu = 0.d0
+  !do i_state = 1, N_states
+  !do i = 1, N_det_ref
+  !  write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state)
+  ! do j = 1, N_det_ref
+  !  accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
+  !  delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+  ! enddo
+  !enddo
+  !second_order_pt_new_2h2p(i_state) = accu(i_state) 
+  !enddo
+  !print*, '2h2p   = ',accu
+   
 
  ! total  
  accu = 0.d0
@@ -234,6 +282,8 @@ END_PROVIDER
  implicit none
  integer :: i,j,i_state
  double precision :: hij
+ double precision :: accu(N_states)
+ accu = 0.d0
  do i_state = 1, N_states
   do i = 1,N_det_ref
    do j = 1,N_det_ref
@@ -241,14 +291,16 @@ END_PROVIDER
     Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) =  hij & 
                                             + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) )
 !   Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) =  Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) 
+    accu(i_State) += psi_ref_coef(i,i_State) * Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) * psi_ref_coef(j,i_State) 
    enddo
   enddo
  enddo
+ print*, 'accu = ',accu + nuclear_repulsion
  END_PROVIDER 
 
-  BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ]
+  BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag_heff) ]
- &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag) ]
+ &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag_heff) ]
- &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ]
+ &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag_heff) ]
   BEGIN_DOC
   ! Eigenvectors/values of the CI matrix
   END_DOC
@@ -269,14 +321,14 @@ END_PROVIDER
   double precision :: overlap(N_det_ref)
   double precision, allocatable :: psi_tmp(:)
   
-  ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors 
+  ! Guess values for the "N_states_diag_heff" states of the CI_dressed_pt2_new_eigenvectors 
   do j=1,min(N_states,N_det_ref)
     do i=1,N_det_ref
       CI_dressed_pt2_new_eigenvectors(i,j) = psi_ref_coef(i,j)
     enddo
   enddo
 
-  do j=min(N_states,N_det_ref)+1,N_states_diag
+  do j=min(N_states,N_det_ref)+1,N_states_diag_heff
     do i=1,N_det_ref
       CI_dressed_pt2_new_eigenvectors(i,j) = 0.d0
     enddo
@@ -408,13 +460,13 @@ END_PROVIDER
           print*,''
           print*,'!!!!!!!!   WARNING  !!!!!!!!!'
           print*,'  Within the ',N_det_ref,'determinants selected'
-          print*,'  and the ',N_states_diag,'states requested'
+          print*,'  and the ',N_states_diag_heff,'states requested'
           print*,'  We did not find any state with S^2 values close to ',expected_s2
           print*,'  We will then set the first N_states eigenvectors of the H matrix'
           print*,'  as the CI_dressed_pt2_new_eigenvectors'
           print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
           print*,''
-          do j=1,min(N_states_diag,N_det_ref)
+          do j=1,min(N_states_diag_heff,N_det_ref)
             do i=1,N_det_ref
               CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
             enddo
@@ -426,8 +478,8 @@ END_PROVIDER
         deallocate(s2_eigvalues)
       else
         call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det_ref,psi_det,N_int,&
-           min(N_det_ref,N_states_diag),size(eigenvectors,1))
+           min(N_det_ref,N_states_diag_heff),size(eigenvectors,1))
-        ! Select the "N_states_diag" states of lowest energy
+        ! Select the "N_states_diag_heff" states of lowest energy
         do j=1,min(N_det_ref,N_states)
           do i=1,N_det_ref
             CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j)
@@ -444,7 +496,7 @@ END_PROVIDER
 END_PROVIDER
  
 
-BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
+BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag_heff) ]
   implicit none
   BEGIN_DOC
   ! N_states lowest eigenvalues of the CI matrix
@@ -453,7 +505,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ]
   integer                        :: j
   character*(8)                  :: st
   call write_time(output_determinants)
-  do j=1,N_states_diag
+  do j=1,N_states_diag_heff
     CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
     write(st,'(I4)') j
     call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))

plugins/MRPT_Utils/new_way.irp.f

@@ -499,9 +499,9 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
       do r = 1, n_virt_orb    ! First virtual
        rorb = list_virt(r) 
        do ispin = 1, 2  ! spin of the couple a-a^dagger (i,r)
-         do state_target = 1, N_states
+        !do state_target = 1, N_states
-          coef_det_pert(i,r,ispin,state_target,1)  += coef_det_pert(i,r,ispin,state_target,2)
+        ! coef_det_pert(i,r,ispin,state_target,1)  += coef_det_pert(i,r,ispin,state_target,2)
-         enddo
+        !enddo
         
         do inint = 1, N_int
          det_tmp(inint,1) = det_pert(inint,1,i,r,ispin)
@@ -509,34 +509,34 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
         enddo
         do jdet = 1, idx(0)
 !      
-         if(idx(jdet).ne.idet)then
-          call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
-          if (exc(0,1,1) == 1) then
-            ! Mono alpha
-            aorb  = (exc(1,2,1))   !!! a^{\dagger}_a 
-            borb  = (exc(1,1,1))   !!! a_{b}
-            jspin = 1
-          else
-            aorb  = (exc(1,2,2))   !!!  a^{\dagger}_a
-            borb  = (exc(1,1,2))   !!!  a_{b}
-            jspin = 2
-          endif
-          
-          call get_excitation_degree(psi_ref(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
-          if(degree_scalar .ne. 2)then
-           print*, 'pb !!!'
-           print*, degree_scalar
-           call debug_det(psi_ref(1,1,idx(jdet)),N_int)
-           call debug_det(det_tmp,N_int)
-           stop
-          endif
-          call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
           double precision :: hij_test
-          hij_test = 0.d0
+         if(idx(jdet).ne.idet)then
-          call  i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test)
+        ! call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int)
-          do state_target = 1, N_states
+        ! if (exc(0,1,1) == 1) then
-           matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
+        !   ! Mono alpha
-          enddo
+        !   aorb  = (exc(1,2,1))   !!! a^{\dagger}_a 
+        !   borb  = (exc(1,1,1))   !!! a_{b}
+        !   jspin = 1
+        ! else
+        !   aorb  = (exc(1,2,2))   !!!  a^{\dagger}_a
+        !   borb  = (exc(1,1,2))   !!!  a_{b}
+        !   jspin = 2
+        ! endif
+        ! 
+        ! call get_excitation_degree(psi_ref(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
+        ! if(degree_scalar .ne. 2)then
+        !  print*, 'pb !!!'
+        !  print*, degree_scalar
+        !  call debug_det(psi_ref(1,1,idx(jdet)),N_int)
+        !  call debug_det(det_tmp,N_int)
+        !  stop
+        ! endif
+        ! call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
+        ! hij_test = 0.d0
+        ! call  i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test)
+        ! do state_target = 1, N_states
+        !  matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
+        ! enddo
          else
           hij_test = 0.d0
           call  i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hij_test)

plugins/MRPT_Utils/psi_active_prov.irp.f

@@ -468,7 +468,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
   endif
  else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then
   do i = 1, N_states
-   delta_e_act(i_state) = -10000000.d0
+   delta_e_act(i_state) = -1.d12
   enddo
  endif
 

plugins/Psiref_CAS/psi_ref.irp.f

@@ -67,3 +67,14 @@ END_PROVIDER
 
 END_PROVIDER
 
+
+ BEGIN_PROVIDER [double precision, electronic_psi_ref_average_value, (N_states)]
+&BEGIN_PROVIDER  [double precision, psi_ref_average_value, (N_states)]
+ implicit none
+ integer :: i,j
+ call u_0_H_u_0(electronic_psi_ref_average_value,psi_ref_coef,N_det_ref,psi_ref,N_int,N_states,psi_det_size) 
+ do i = 1, N_states
+  psi_ref_average_value(i) = electronic_psi_ref_average_value(i) + nuclear_repulsion
+ enddo
+
+END_PROVIDER

plugins/loc_cele/loc.f

@@ -18,7 +18,7 @@ C
 
       zprt=.true.
       niter=1000000
-      conv=1.d-8
+      conv=1.d-10
 
 C      niter=1000000
 C      conv=1.d-6

plugins/loc_cele/loc_cele.irp.f

@@ -101,10 +101,27 @@
        cmoref = 0.d0
        irot = 0
 
-       irot(1,1) = 11
+       irot(1,1) = 5 
-       irot(2,1) = 12
+       irot(2,1) = 6 
-       cmoref(15,1,1)   = 1.d0   ! 
+       cmoref(6,1,1)   = 1d0 
-       cmoref(14,2,1)   = 1.d0   ! 
+       cmoref(26,2,1)   = 1d0 
+
+!       !!! H2O
+!      irot(1,1) = 4 
+!      irot(2,1) = 5 
+!      irot(3,1) = 6 
+!      irot(4,1) = 7 
+!      ! O pz 
+!      cmoref(5,1,1)   = 1.55362d0 
+!      cmoref(6,1,1)   = 1.07578d0 
+
+!      cmoref(5,2,1)   =  1.55362d0 
+!      cmoref(6,2,1)   = -1.07578d0 
+!      ! O px - pz 
+!      ! H1 
+!      cmoref(16,3,1)   = 1.d0
+!      ! H1 
+!      cmoref(21,4,1)   = 1.d0
 
 ! ESATRIENE with 3 bonding and anti bonding orbitals 
 ! First bonding orbital for esa
@@ -150,19 +167,19 @@
 ! ESATRIENE with 1 central bonding and anti bonding orbitals 
 ! AND 4 radical orbitals
 ! First radical orbital 
-       cmoref(7,1,1)   = 1.d0   ! 
+!      cmoref(7,1,1)   = 1.d0   ! 
 ! Second radical orbital 
-       cmoref(26,2,1)  = 1.d0   ! 
+!      cmoref(26,2,1)  = 1.d0   ! 
 ! First bonding orbital 
-       cmoref(45,3,1)  = 1.d0   ! 
+!      cmoref(45,3,1)  = 1.d0   ! 
-       cmoref(64,3,1)  = 1.d0   ! 
+!      cmoref(64,3,1)  = 1.d0   ! 
 ! Third radical orbital for esa
-       cmoref(83,4,1)  = 1.d0   ! 
+!      cmoref(83,4,1)  = 1.d0   ! 
 ! Fourth radical orbital for esa
-       cmoref(102,5,1) = 1.d0   ! 
+!      cmoref(102,5,1) = 1.d0   ! 
 ! First anti bonding orbital 
-       cmoref(45,6,1)  = 1.d0   ! 
+!      cmoref(45,6,1)  = 1.d0   ! 
-       cmoref(64,6,1)  =-1.d0   ! 
+!      cmoref(64,6,1)  =-1.d0   ! 
 
 
        do i = 1, nrot(1)

plugins/loc_cele/loc_exchange_int_act.irp.f

@@ -19,16 +19,17 @@ program loc_int
   do j = i+1, n_act_orb
    jorb = list_act(j)
    iorder(jorb) = jorb
-   exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   if(list_core_inact_check(jorb) == .False.)then
+    exchange_int(jorb) = 0.d0
+   else 
+    exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb)
+   endif
   enddo
   n_rot += 1
   call dsort(exchange_int,iorder,mo_tot_num)
   indices(n_rot,1) = iorb
   indices(n_rot,2) = iorder(1)
   list_core_inact_check(iorder(1)) = .False.
-  print*,indices(n_rot,1),indices(n_rot,2)
-  print*,''
-  print*,''
  enddo
   print*,'****************************'
   print*,'-+++++++++++++++++++++++++'

src/Integrals_Monoelec/pot_ao_ints.irp.f

@@ -185,7 +185,7 @@ include 'Utils/constants.include.F'
   enddo
   const_factor = dist*rho
   const = p * dist_integral
-  if(const_factor > 80.d0)then
+  if(const_factor > 1000.d0)then
    NAI_pol_mult = 0.d0
    return
   endif