diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f index e5d925a3..af3713c5 100644 --- a/plugins/MRPT/MRPT_Utils.main.irp.f +++ b/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+PT2 ', i,' = ', CI_energy(i)+second_order_pt_new(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,' = ', 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 diff --git a/plugins/MRPT/print_1h2p.irp.f b/plugins/MRPT/print_1h2p.irp.f index 2739340b..85ddcda8 100644 --- a/plugins/MRPT/print_1h2p.irp.f +++ b/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_act_orb - write(*,'(I3,x,100(F16.10,X))')i,one_creat(i,:,1) -! write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(1,4,1,2,1) -! + do i =1, n_core_inact_orb + write(*,'(I3,x,100(F16.10,X))')list_core_inact(i),fock_core_inactive_total_spin_trace(list_core_inact(i),1) + enddo + 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 diff --git a/plugins/MRPT_Utils/EZFIO.cfg b/plugins/MRPT_Utils/EZFIO.cfg index 948aa735..db7b127a 100644 --- a/plugins/MRPT_Utils/EZFIO.cfg +++ b/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 + diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index 8f6a7eb6..f8782bec 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/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 diff --git a/plugins/MRPT_Utils/ezfio_interface.irp.f b/plugins/MRPT_Utils/ezfio_interface.irp.f index ebe0bf52..91c7ea63 100644 --- a/plugins/MRPT_Utils/ezfio_interface.irp.f +++ b/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 diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index f241d35a..1fd8cb03 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/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 - !do i_state = 1, N_states - ! delta_e(i_state) = 1.d+20 - !enddo + + ! !!!!!!!!!!!!! SHIFTED BK + ! double precision :: hjj + ! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj) + ! delta_e(1) = CI_electronic_energy(1) - hjj + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + endif hij_array(index_i) = hialpha do i_state = 1,N_states diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 09efc536..35940404 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/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 @@ -167,6 +199,22 @@ second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) 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 @@ -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_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (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_heff) ] + &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)) - ! Select the "N_states_diag" states of lowest energy + min(N_det_ref,N_states_diag_heff),size(eigenvectors,1)) + ! 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)) diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index dc921551..3cfa7154 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/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 - coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2) - enddo + !do state_target = 1, N_states + ! coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2) + !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 - 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 + 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) + ! 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) diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index bd31dc1d..95f7479e 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/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 diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index d3b6c28f..0729a540 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/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 diff --git a/plugins/loc_cele/loc.f b/plugins/loc_cele/loc.f index edc3aa7a..ed8b9a76 100644 --- a/plugins/loc_cele/loc.f +++ b/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 diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f index 2d47c633..2dda522e 100644 --- a/plugins/loc_cele/loc_cele.irp.f +++ b/plugins/loc_cele/loc_cele.irp.f @@ -101,10 +101,27 @@ cmoref = 0.d0 irot = 0 - irot(1,1) = 11 - irot(2,1) = 12 - cmoref(15,1,1) = 1.d0 ! - cmoref(14,2,1) = 1.d0 ! + irot(1,1) = 5 + irot(2,1) = 6 + cmoref(6,1,1) = 1d0 + 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(64,3,1) = 1.d0 ! +! cmoref(45,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(64,6,1) =-1.d0 ! +! cmoref(45,6,1) = 1.d0 ! +! cmoref(64,6,1) =-1.d0 ! do i = 1, nrot(1) diff --git a/plugins/loc_cele/loc_exchange_int_act.irp.f b/plugins/loc_cele/loc_exchange_int_act.irp.f index f332dd5d..c4dcf75c 100644 --- a/plugins/loc_cele/loc_exchange_int_act.irp.f +++ b/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*,'-+++++++++++++++++++++++++' diff --git a/src/Integrals_Monoelec/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f index 7116d2c7..aef8a060 100644 --- a/src/Integrals_Monoelec/pot_ao_ints.irp.f +++ b/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