diff --git a/README.md b/README.md index e313f444..5372b7ac 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,7 @@ -Quantum package -=============== - +![QP](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/qp.png) [![Build Status](https://travis-ci.org/LCPQ/quantum_package.svg?branch=master)](https://travis-ci.org/LCPQ/quantum_package) - [![Gitter](https://badges.gitter.im/Join%20Chat.svg)](https://gitter.im/LCPQ/quantum_package?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) - Set of quantum chemistry programs and libraries. (under GNU GENERAL PUBLIC LICENSE v2) diff --git a/config/ifort.cfg b/config/ifort.cfg index 2b2fe0a2..cb6dc1ef 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -31,14 +31,15 @@ OPENMP : 1 ; Append OpenMP flags # -ftz : Flushes denormal results to zero # [OPT] -FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g +FC : -traceback +FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g -traceback # Profiling flags ################# # [PROFILE] -FC : -p -g -FCFLAGS : -xSSE4.2 -O2 -ip -ftz +FC : -p -g -traceback +FCFLAGS : -xSSE4.2 -O2 -ip -ftz # Debugging flags ################# @@ -51,13 +52,13 @@ FCFLAGS : -xSSE4.2 -O2 -ip -ftz # [DEBUG] FC : -g -traceback -FCFLAGS : -xSSE2 -C -fpe0 +FCFLAGS : -xSSE2 -C IRPF90_FLAGS : --openmp # OpenMP flags ################# # [OPENMP] -FC : -openmp +FC : -qopenmp IRPF90_FLAGS : --openmp diff --git a/data/qp.png b/data/qp.png new file mode 100644 index 00000000..777e5ac0 Binary files /dev/null and b/data/qp.png differ diff --git a/plugins/All_singles/H_apply.irp.f b/plugins/All_singles/H_apply.irp.f index d0a41f90..f34f003c 100644 --- a/plugins/All_singles/H_apply.irp.f +++ b/plugins/All_singles/H_apply.irp.f @@ -8,10 +8,9 @@ s.unset_skip() s.filter_only_1h1p() print s -s = H_apply("just_mono") +s = H_apply("just_mono",do_double_exc=False) s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() -s.unset_double_excitations() print s END_SHELL diff --git a/plugins/All_singles/all_singles.irp.f b/plugins/All_singles/all_singles.irp.f index 3b5c5cce..ad8648c7 100644 --- a/plugins/All_singles/all_singles.irp.f +++ b/plugins/All_singles/all_singles.irp.f @@ -15,7 +15,7 @@ subroutine routine integer :: N_st, degree double precision,allocatable :: E_before(:) integer :: n_det_before - N_st = N_states + N_st = N_states_diag allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) i = 0 print*,'N_det = ',N_det diff --git a/plugins/CAS_SD/H_apply.irp.f b/plugins/CAS_SD/H_apply.irp.f index aa393bc7..35c45fb6 100644 --- a/plugins/CAS_SD/H_apply.irp.f +++ b/plugins/CAS_SD/H_apply.irp.f @@ -20,22 +20,18 @@ print s s = H_apply("CAS_S",do_double_exc=False) -s.unset_double_excitations() print s s = H_apply("CAS_S_selected_no_skip",do_double_exc=False) -s.unset_double_excitations() s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() print s s = H_apply("CAS_S_selected",do_double_exc=False) -s.unset_double_excitations() s.set_selection_pt2("epstein_nesbet_2x2") print s s = H_apply("CAS_S_PT2",do_double_exc=False) -s.unset_double_excitations() s.set_perturbation("epstein_nesbet_2x2") print s diff --git a/plugins/DDCI_selected/ddci.irp.f b/plugins/DDCI_selected/ddci.irp.f index 248671b1..0bfb324f 100644 --- a/plugins/DDCI_selected/ddci.irp.f +++ b/plugins/DDCI_selected/ddci.irp.f @@ -3,10 +3,10 @@ program ddci integer :: i,k - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) integer :: N_st, degree - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + N_st = N_states_diag + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) character*(64) :: perturbation pt2 = 1.d0 @@ -27,6 +27,8 @@ program ddci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_DDCI_selection(pt2, norm_pert, H_pert_diag, N_st) @@ -47,8 +49,21 @@ program ddci print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 + print *, 'E+PT2 = ', E_before+pt2 print *, '-----' + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy call ezfio_set_ddci_selected_energy(CI_energy) enddo if(do_pt2_end)then diff --git a/plugins/FOBOCI/EZFIO.cfg b/plugins/FOBOCI/EZFIO.cfg index d4a10add..88189608 100644 --- a/plugins/FOBOCI/EZFIO.cfg +++ b/plugins/FOBOCI/EZFIO.cfg @@ -1,6 +1,13 @@ -[threshold_singles] +[threshold_lmct] type: double precision -doc: threshold to select the pertinent single excitations at second order +doc: threshold to select the pertinent LMCT excitations at second order +interface: ezfio,provider,ocaml +default: 0.01 + + +[threshold_mlct] +type: double precision +doc: threshold to select the pertinent MLCT excitations at second order interface: ezfio,provider,ocaml default: 0.01 @@ -16,6 +23,20 @@ doc: if true, you do the FOBOCI calculation perturbatively interface: ezfio,provider,ocaml default: .False. + +[speed_up_convergence_foboscf] +type: logical +doc: if true, the threshold of the FOBO-SCF algorithms are increased with the iterations +interface: ezfio,provider,ocaml +default: .True. + + +[dressing_2h2p] +type: logical +doc: if true, you do dress with 2h2p excitations each FOBOCI matrix +interface: ezfio,provider,ocaml +default: .False. + [second_order_h] type: logical doc: if true, you do the FOBOCI calculation using second order intermediate Hamiltonian diff --git a/plugins/FOBOCI/H_apply.irp.f b/plugins/FOBOCI/H_apply.irp.f index 0a488753..d8ab02f1 100644 --- a/plugins/FOBOCI/H_apply.irp.f +++ b/plugins/FOBOCI/H_apply.irp.f @@ -18,8 +18,22 @@ print s -s = H_apply("standard") +s = H_apply("only_1h2p") s.set_selection_pt2("epstein_nesbet") +s.filter_only_1h2p() +s.unset_skip() +print s + +s = H_apply("only_2h2p") +s.set_selection_pt2("epstein_nesbet") +s.filter_only_2h2p() +s.unset_skip() +print s + + +s = H_apply("only_2p") +s.set_selection_pt2("epstein_nesbet") +s.filter_only_2p() s.unset_skip() print s diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES index adeefe99..f6c0c1c4 100644 --- a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES +++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Generators_restart Selectors_no_sorted +Perturbation Selectors_no_sorted Hartree_Fock diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f index e2c4c01e..0594e56e 100644 --- a/plugins/FOBOCI/all_singles.irp.f +++ b/plugins/FOBOCI/all_singles.irp.f @@ -6,9 +6,9 @@ subroutine all_single double precision,allocatable :: E_before(:) N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) - selection_criterion = 1.d-8 + selection_criterion = 0.d0 soft_touch selection_criterion - threshold_davidson = 1.d-5 + threshold_davidson = 1.d-9 soft_touch threshold_davidson davidson_criterion i = 0 print*,'Doing all the mono excitations !' @@ -52,10 +52,173 @@ subroutine all_single enddo endif E_before = CI_energy + !!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO + exit enddo - threshold_davidson = 1.d-10 +! threshold_davidson = 1.d-8 +! soft_touch threshold_davidson davidson_criterion +! call diagonalize_CI + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + do i = 1, max(2,N_det_generators) + print*,'psi_coef = ',psi_coef(i,1) + enddo + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_1h2p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion + threshold_davidson = 1.d-5 soft_touch threshold_davidson davidson_criterion - call diagonalize_CI + i = 0 + print*,'' + print*,'' + print*,'' + print*,'' + print*,'' + print*,'*****************************' + print*,'Doing all the 1h2P excitations' + print*,'*****************************' + print*,'' + print*,'' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + i = 0 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_only_1h2p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + + enddo + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_2h2p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'' + print*,'' + print*,'' + print*,'' + print*,'' + print*,'*****************************' + print*,'Doing all the 2h2P excitations' + print*,'*****************************' + print*,'' + print*,'' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + i = 0 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_only_2h2p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + + enddo print*,'Final Step ' print*,'N_det = ',N_det do i = 1, N_states_diag @@ -67,10 +230,89 @@ subroutine all_single do i = 1, 2 print*,'psi_coef = ',psi_coef(i,1) enddo -! call save_wavefunction deallocate(pt2,norm_pert,E_before) end +subroutine all_2p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'' + print*,'' + print*,'' + print*,'' + print*,'' + print*,'*****************************' + print*,'Doing all the 2P excitations' + print*,'*****************************' + print*,'' + print*,'' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + i = 0 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_only_2p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + + enddo + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + deallocate(pt2,norm_pert,E_before) + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo +end + subroutine all_single_no_1h_or_1p implicit none integer :: i,k @@ -79,6 +321,8 @@ subroutine all_single_no_1h_or_1p double precision,allocatable :: E_before(:) N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion threshold_davidson = 1.d-5 soft_touch threshold_davidson davidson_criterion i = 0 @@ -124,7 +368,7 @@ subroutine all_single_no_1h_or_1p endif E_before = CI_energy enddo - threshold_davidson = 1.d-10 + threshold_davidson = 1.d-16 soft_touch threshold_davidson davidson_criterion call diagonalize_CI print*,'Final Step ' @@ -215,85 +459,6 @@ subroutine all_single_no_1h_or_1p_or_2p deallocate(pt2,norm_pert,E_before) end - -subroutine all_2p - implicit none - integer :: i,k - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer :: N_st, degree - double precision,allocatable :: E_before(:) - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) - selection_criterion = 0.d0 - soft_touch selection_criterion - threshold_davidson = 1.d-5 - soft_touch threshold_davidson davidson_criterion - i = 0 - print*,'' - print*,'' - print*,'' - print*,'' - print*,'' - print*,'*****************************' - print*,'Doing all the 2P excitations' - print*,'*****************************' - print*,'' - print*,'' - print*,'N_det = ',N_det - print*,'n_det_max = ',n_det_max - print*,'pt2_max = ',pt2_max - print*,'N_det_generators = ',N_det_generators - pt2=-1.d0 - E_before = ref_bitmask_energy - - print*,'Initial Step ' - print*,'Inital determinants ' - print*,'N_det = ',N_det - do i = 1, N_states_diag - print*,'' - print*,'i = ',i - print*,'E = ',CI_energy(i) - print*,'S^2 = ',CI_eigenvectors_s2(i) - enddo - n_det_max = 100000 - i = 0 - do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) - i += 1 - print*,'-----------------------' - print*,'i = ',i - call H_apply_standard(pt2, norm_pert, H_pert_diag, N_st) - call diagonalize_CI - print*,'N_det = ',N_det - print*,'E = ',CI_energy(1) - print*,'pt2 = ',pt2(1) - print*,'E+PT2 = ',E_before + pt2(1) - if(N_states_diag.gt.1)then - print*,'Variational Energy difference' - do i = 2, N_st - print*,'Delta E = ',CI_energy(i) - CI_energy(1) - enddo - endif - if(N_states.gt.1)then - print*,'Variational + perturbative Energy difference' - do i = 2, N_st - print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) - enddo - endif - E_before = CI_energy - - enddo - print*,'Final Step ' - print*,'N_det = ',N_det - do i = 1, N_states_diag - print*,'' - print*,'i = ',i - print*,'E = ',CI_energy(i) - print*,'S^2 = ',CI_eigenvectors_s2(i) - enddo -! call save_wavefunction - deallocate(pt2,norm_pert,E_before) -end - subroutine all_1h_1p_routine implicit none integer :: i,k diff --git a/plugins/FOBOCI/all_singles_split.irp.f b/plugins/FOBOCI/all_singles_split.irp.f index e7b0943f..9ddf369a 100644 --- a/plugins/FOBOCI/all_singles_split.irp.f +++ b/plugins/FOBOCI/all_singles_split.irp.f @@ -5,7 +5,7 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators_input) double precision, intent(inout) :: dressing_matrix(Ndet_generators_input,Ndet_generators_input) double precision, intent(in) :: psi_coef_generators_input(ndet_generators_input,n_states) - integer :: i,i_hole + integer :: i,i_hole,j n_det_max_jacobi = 50 soft_touch n_det_max_jacobi do i = 1, n_inact_orb @@ -22,56 +22,339 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) call all_single - threshold_davidson = 1.d-10 - soft_touch threshold_davidson davidson_criterion - call diagonalize_CI +! call diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) enddo + + do i = 1, n_act_orb + i_hole = list_act(i) + print*,'' + print*,'Doing all the single excitations from the orbital ' + print*,i_hole + print*,'' + print*,'' + threshold_davidson = 1.d-4 + soft_touch threshold_davidson davidson_criterion + call modify_bitmasks_for_hole(i_hole) + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call all_single +! call diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) + call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) + enddo + + do i = 1, n_virt_orb + i_hole = list_virt(i) + print*,'' + print*,'Doing all the single excitations from the orbital ' + print*,i_hole + print*,'' + print*,'' + threshold_davidson = 1.d-4 + soft_touch threshold_davidson davidson_criterion + call modify_bitmasks_for_hole(i_hole) + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call all_single +! call diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) + call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) + enddo + n_det_max_jacobi = 1000 soft_touch n_det_max_jacobi end -subroutine all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + +subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) implicit none use bitmasks + integer, intent(in) :: i_particl double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) - double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) - integer :: i,i_hole + double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) + integer :: i,j n_det_max_jacobi = 50 soft_touch n_det_max_jacobi - integer :: n_det_1h1p,n_det_2h1p - integer(bit_kind), allocatable :: psi_ref_out(:,:,:) - integer(bit_kind), allocatable :: psi_1h1p(:,:,:) - integer(bit_kind), allocatable :: psi_2h1p(:,:,:) - double precision, allocatable :: psi_ref_coef_out(:,:) - double precision, allocatable :: psi_coef_1h1p(:,:) - double precision, allocatable :: psi_coef_2h1p(:,:) - call all_single_no_1h_or_1p + call all_single threshold_davidson = 1.d-12 soft_touch threshold_davidson davidson_criterion call diagonalize_CI - call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) - allocate(psi_ref_out(N_int,2,N_det_generators)) - allocate(psi_1h1p(N_int,2,n_det_1h1p)) - allocate(psi_2h1p(N_int,2,n_det_2h1p)) - allocate(psi_ref_coef_out(N_det_generators,N_states)) - allocate(psi_coef_1h1p(n_det_1h1p,N_states)) - allocate(psi_coef_2h1p(n_det_2h1p,N_states)) - call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) - call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_1h1p,psi_coef_1h1p,n_det_1h1p) - call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_2h1p,psi_coef_2h1p,n_det_2h1p) - deallocate(psi_ref_out) - deallocate(psi_1h1p) - deallocate(psi_2h1p) - deallocate(psi_ref_coef_out) - deallocate(psi_coef_1h1p) - deallocate(psi_coef_2h1p) + + + double precision, allocatable :: matrix_ref_1h_1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h2p(:,:) + double precision, allocatable :: psi_coef_ref_1h_1p(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_1h2p(:,:) + integer(bit_kind), allocatable :: psi_det_1h2p(:,:,:) + integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:) + integer :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + double precision :: hka + double precision,allocatable :: eigenvectors(:,:), eigenvalues(:) + + + call give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p) + + allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h2p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)) + allocate(psi_det_1h2p(N_int,2,n_det_1h2p), psi_coef_1h2p(n_det_1h2p,N_states)) + allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states)) + + call give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p) + + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka) + matrix_ref_1h_1p(i,j) = hka + enddo + enddo + matrix_ref_1h_1p_dressing_1h1p = 0.d0 + matrix_ref_1h_1p_dressing_1h2p = 0.d0 + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h2p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h2p,psi_coef_1h2p,n_det_1h2p) + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h1p,psi_coef_1h1p,n_det_1h1p) + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_1h2p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j) + enddo + enddo + + allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p)) + call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p) +!do j = 1, n_det_ref_1h_1p +! print*,'coef = ',eigenvectors(j,1) +!enddo + print*,'' + print*,'-----------------------' + print*,'-----------------------' + print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion + print*,'-----------------------' + ! Extract the + integer, allocatable :: index_generator(:) + integer :: n_det_generators_tmp,degree + n_det_generators_tmp = 0 + allocate(index_generator(n_det_ref_1h_1p)) + do i = 1, n_det_ref_1h_1p + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(1,1,i), degree, N_int) + if(degree == 0)then + n_det_generators_tmp +=1 + index_generator(n_det_generators_tmp) = i + endif + enddo + enddo + if(n_det_generators_tmp .ne. n_det_generators)then + print*,'PB !!!' + print*,'if(n_det_generators_tmp .ne. n_det_genrators)then' + stop + endif + do i = 1, N_det_generators + print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1) + do j = 1, N_det_generators + dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j)) + dressing_matrix_1h2p(i,j) += matrix_ref_1h_1p_dressing_1h2p(index_generator(i),index_generator(j)) + enddo + enddo + print*,'-----------------------' + print*,'-----------------------' + + + deallocate(matrix_ref_1h_1p) + deallocate(matrix_ref_1h_1p_dressing_1h1p) + deallocate(matrix_ref_1h_1p_dressing_1h2p) + deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p) + deallocate(psi_det_1h2p, psi_coef_1h2p) + deallocate(psi_det_1h1p, psi_coef_1h1p) + deallocate(eigenvectors,eigenvalues) + deallocate(index_generator) + + +end + +subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p) + implicit none + use bitmasks + integer, intent(in) :: i_hole + double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) + integer :: i,j + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + call all_single + + threshold_davidson = 1.d-12 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + + + + double precision, allocatable :: matrix_ref_1h_1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_2h1p(:,:) + double precision, allocatable :: psi_coef_ref_1h_1p(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_2h1p(:,:) + integer(bit_kind), allocatable :: psi_det_2h1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:) + integer :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p + double precision :: hka + double precision,allocatable :: eigenvectors(:,:), eigenvalues(:) + + + call give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p) + + allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_2h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)) + allocate(psi_det_2h1p(N_int,2,n_det_2h1p), psi_coef_2h1p(n_det_2h1p,N_states)) + allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states)) + + call give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p) + + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka) + matrix_ref_1h_1p(i,j) = hka + enddo + enddo + matrix_ref_1h_1p_dressing_1h1p = 0.d0 + matrix_ref_1h_1p_dressing_2h1p = 0.d0 + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_2h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_2h1p,psi_coef_2h1p,n_det_2h1p) + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h1p,psi_coef_1h1p,n_det_1h1p) + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j) + enddo + enddo + + allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p)) + call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p) +!do j = 1, n_det_ref_1h_1p +! print*,'coef = ',eigenvectors(j,1) +!enddo + print*,'' + print*,'-----------------------' + print*,'-----------------------' + print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion + print*,'-----------------------' + ! Extract the + integer, allocatable :: index_generator(:) + integer :: n_det_generators_tmp,degree + n_det_generators_tmp = 0 + allocate(index_generator(n_det_ref_1h_1p)) + do i = 1, n_det_ref_1h_1p + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(1,1,i), degree, N_int) + if(degree == 0)then + n_det_generators_tmp +=1 + index_generator(n_det_generators_tmp) = i + endif + enddo + enddo + if(n_det_generators_tmp .ne. n_det_generators)then + print*,'PB !!!' + print*,'if(n_det_generators_tmp .ne. n_det_genrators)then' + stop + endif + do i = 1, N_det_generators + print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1) + do j = 1, N_det_generators + dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j)) + dressing_matrix_2h1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(index_generator(i),index_generator(j)) + enddo + enddo + print*,'-----------------------' + print*,'-----------------------' + + + deallocate(matrix_ref_1h_1p) + deallocate(matrix_ref_1h_1p_dressing_1h1p) + deallocate(matrix_ref_1h_1p_dressing_2h1p) + deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p) + deallocate(psi_det_2h1p, psi_coef_2h1p) + deallocate(psi_det_1h1p, psi_coef_1h1p) + deallocate(eigenvectors,eigenvalues) + deallocate(index_generator) +!return +! + +!integer(bit_kind), allocatable :: psi_ref_out(:,:,:) +!integer(bit_kind), allocatable :: psi_1h1p(:,:,:) +!integer(bit_kind), allocatable :: psi_2h1p(:,:,:) +!integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) +!double precision, allocatable :: psi_ref_coef_out(:,:) +!double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) + +!call all_single_no_1h_or_1p + +!call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p) +!allocate(psi_ref_out(N_int,2,N_det_generators)) +!allocate(psi_1h1p(N_int,2,n_det_1h1p)) +!allocate(psi_2h1p(N_int,2,n_det_2h1p)) +!allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) +!allocate(psi_ref_coef_out(N_det_generators,N_states)) +!allocate(psi_coef_1h1p(n_det_1h1p,N_states)) +!allocate(psi_coef_2h1p(n_det_2h1p,N_states)) +!allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) +!call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) +!do i = 1, n_det_extra_1h_or_1p +! print*,'----' +! print*,'c = ',psi_coef_extra_1h_or_1p(i,1) +! call debug_det(psi_extra_1h_or_1p(1,1,i),N_int) +! print*,'----' +!enddo +!call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_1h1p,psi_coef_1h1p,n_det_1h1p) +!print*,'Dressing 1h1p ' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_1h1p(j,:) +!enddo + +!call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_2h1p,psi_coef_2h1p,n_det_2h1p) +!print*,'Dressing 2h1p ' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_2h1p(j,:) +!enddo + +!call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) +!print*,',dressing_matrix_extra_1h_or_1p' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_extra_1h_or_1p(j,:) +!enddo + + +!deallocate(psi_ref_out) +!deallocate(psi_1h1p) +!deallocate(psi_2h1p) +!deallocate(psi_extra_1h_or_1p) +!deallocate(psi_ref_coef_out) +!deallocate(psi_coef_1h1p) +!deallocate(psi_coef_2h1p) +!deallocate(psi_coef_extra_1h_or_1p) end @@ -197,47 +480,56 @@ subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) soft_touch n_det_max_jacobi end -subroutine all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) - implicit none - use bitmasks - double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) - double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) - integer :: i,i_hole - n_det_max_jacobi = 50 - soft_touch n_det_max_jacobi - - integer :: n_det_1h1p,n_det_1h2p - integer(bit_kind), allocatable :: psi_ref_out(:,:,:) - integer(bit_kind), allocatable :: psi_1h1p(:,:,:) - integer(bit_kind), allocatable :: psi_1h2p(:,:,:) - double precision, allocatable :: psi_ref_coef_out(:,:) - double precision, allocatable :: psi_coef_1h1p(:,:) - double precision, allocatable :: psi_coef_1h2p(:,:) - call all_single_no_1h_or_1p_or_2p - - threshold_davidson = 1.d-12 - soft_touch threshold_davidson davidson_criterion - call diagonalize_CI - call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) - allocate(psi_ref_out(N_int,2,N_det_generators)) - allocate(psi_1h1p(N_int,2,n_det_1h1p)) - allocate(psi_1h2p(N_int,2,n_det_1h2p)) - allocate(psi_ref_coef_out(N_det_generators,N_states)) - allocate(psi_coef_1h1p(n_det_1h1p,N_states)) - allocate(psi_coef_1h2p(n_det_1h2p,N_states)) - call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) - call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_1h1p,psi_coef_1h1p,n_det_1h1p) - call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_1h2p,psi_coef_1h2p,n_det_1h2p) - - deallocate(psi_ref_out) - deallocate(psi_1h1p) - deallocate(psi_1h2p) - deallocate(psi_ref_coef_out) - deallocate(psi_coef_1h1p) - deallocate(psi_coef_1h2p) - -end +! subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) +! implicit none +! use bitmasks +! integer, intent(in ) :: i_particl +! double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) +! double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) +! double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) +! integer :: i +! n_det_max_jacobi = 50 +! soft_touch n_det_max_jacobi +! +! integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p +! integer(bit_kind), allocatable :: psi_ref_out(:,:,:) +! integer(bit_kind), allocatable :: psi_1h1p(:,:,:) +! integer(bit_kind), allocatable :: psi_1h2p(:,:,:) +! integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) +! double precision, allocatable :: psi_ref_coef_out(:,:) +! double precision, allocatable :: psi_coef_1h1p(:,:) +! double precision, allocatable :: psi_coef_1h2p(:,:) +! double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) +!!!!call all_single_no_1h_or_1p_or_2p +! call all_single +! +! threshold_davidson = 1.d-12 +! soft_touch threshold_davidson davidson_criterion +! call diagonalize_CI +! call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p) +! allocate(psi_ref_out(N_int,2,N_det_generators)) +! allocate(psi_1h1p(N_int,2,n_det_1h1p)) +! allocate(psi_1h2p(N_int,2,n_det_1h2p)) +! allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) +! allocate(psi_ref_coef_out(N_det_generators,N_states)) +! allocate(psi_coef_1h1p(n_det_1h1p,N_states)) +! allocate(psi_coef_1h2p(n_det_1h2p,N_states)) +! allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) +! call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) +! call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_1h1p,psi_coef_1h1p,n_det_1h1p) +! call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_1h2p,psi_coef_1h2p,n_det_1h2p) +! call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) +! +! deallocate(psi_ref_out) +! deallocate(psi_1h1p) +! deallocate(psi_1h2p) +! deallocate(psi_ref_coef_out) +! deallocate(psi_coef_1h1p) +! deallocate(psi_coef_1h2p) +! +! end diff --git a/plugins/FOBOCI/collect_all_lmct.irp.f b/plugins/FOBOCI/collect_all_lmct.irp.f new file mode 100644 index 00000000..96eb2858 --- /dev/null +++ b/plugins/FOBOCI/collect_all_lmct.irp.f @@ -0,0 +1,436 @@ +use bitmasks + +subroutine collect_lmct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + print*,'COLLECTING THE PERTINENT LMCT (1h)' + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-2)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_mlct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + print*,'COLLECTING THE PERTINENT MLCT (1p)' + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-3)then + n_couples +=1 + hole_particle(n_couples,1) = iorb + hole_particle(n_couples,2) = jorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_lmct_mlct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + double precision, allocatable :: tmp(:,:) + print*,'COLLECTING THE PERTINENT LMCT (1h)' + print*,'AND THE PERTINENT MLCT (1p)' + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then + n_couples +=1 + hole_particle(n_couples,1) = iorb + hole_particle(n_couples,2) = jorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_1h1p(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + double precision, allocatable :: tmp(:,:) + print*,'COLLECTING THE PERTINENT 1h1p' + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_virt_orb + iorb = list_virt(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-2)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + + +subroutine set_lmct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + print*,'i_hole,i_particle 2 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + print*,'i_hole,i_particle 1 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_mlct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_mlct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + print*,'i_hole,i_particle 2 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + print*,'i_hole,i_particle 1 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_lmct_mlct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct_mlct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_lmct_mlct_to_psi_det + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct_mlct(hole_particle,n_couples) + call set_psi_det_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_generators_restart + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + + N_det = N_det_total + integer :: k + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total)) + enddo + enddo + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) +end + +subroutine set_1h1p_to_psi_det + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_1h1p(hole_particle,n_couples) + call set_psi_det_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_generators_restart + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + + N_det = N_det_total + integer :: k + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total)) + enddo + enddo + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) +end + diff --git a/plugins/FOBOCI/corr_energy_2h2p.irp.f b/plugins/FOBOCI/corr_energy_2h2p.irp.f new file mode 100644 index 00000000..ada46bf2 --- /dev/null +++ b/plugins/FOBOCI/corr_energy_2h2p.irp.f @@ -0,0 +1,425 @@ + BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_ab_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_bb_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_a, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_b, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_double, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_2h2p] + use bitmasks + print*,'' + print*,'Providing the 2h2p correlation energy' + print*,'' + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_2h2p = 0.d0 + corr_energy_2h2p_ab_2_orb = 0.d0 + corr_energy_2h2p_bb_2_orb = 0.d0 + corr_energy_2h2p_per_orb_ab = 0.d0 + corr_energy_2h2p_per_orb_aa = 0.d0 + corr_energy_2h2p_per_orb_bb = 0.d0 + corr_energy_2h2p_for_1h1p_a = 0.d0 + corr_energy_2h2p_for_1h1p_b = 0.d0 + corr_energy_2h2p_for_1h1p_double = 0.d0 + do i = 1, n_inact_orb ! beta + i_hole = list_inact(i) + do k = 1, n_virt_orb ! beta + k_part = list_virt(k) + do j = 1, n_inact_orb ! alpha + j_hole = list_inact(j) + do l = 1, n_virt_orb ! alpha + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = (ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = hij*hij/delta_e + total_corr_e_2h2p += contrib + ! Single orbital contribution + corr_energy_2h2p_per_orb_ab(i_hole) += contrib + corr_energy_2h2p_per_orb_ab(k_part) += contrib + ! Couple of orbital contribution for the single 1h1p + corr_energy_2h2p_for_1h1p_a(j_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_a(l_part,j_hole) += contrib + corr_energy_2h2p_for_1h1p_b(j_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_b(l_part,j_hole) += contrib + ! Couple of orbital contribution for the double 1h1p + corr_energy_2h2p_for_1h1p_double(i_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_double(l_part,i_hole) += contrib + + corr_energy_2h2p_ab_2_orb(i_hole,j_hole) += contrib + corr_energy_2h2p_ab_2_orb(j_hole,i_hole) += contrib + corr_energy_2h2p_ab_2_orb(i_hole,k_part) += contrib + corr_energy_2h2p_ab_2_orb(k_part,i_hole) += contrib + corr_energy_2h2p_ab_2_orb(k_part,l_part) += contrib + corr_energy_2h2p_ab_2_orb(l_part,k_part) += contrib + enddo + enddo + enddo + enddo + + ! alpha alpha correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_virt_orb + k_part = list_virt(k) + do l = k+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h2p += contrib + ! Single orbital contribution + corr_energy_2h2p_per_orb_aa(i_hole) += contrib + corr_energy_2h2p_per_orb_aa(k_part) += contrib + ! Couple of orbital contribution for the single 1h1p + corr_energy_2h2p_for_1h1p_a(i_hole,k_part) += contrib + corr_energy_2h2p_for_1h1p_a(k_part,i_hole) += contrib + enddo + enddo + enddo + enddo + + ! beta beta correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_virt_orb + k_part = list_virt(k) + do l = k+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h2p += contrib + ! Single orbital contribution + corr_energy_2h2p_per_orb_bb(i_hole) += contrib + corr_energy_2h2p_per_orb_bb(k_part) += contrib + corr_energy_2h2p_for_1h1p_b(i_hole,k_part) += contrib + corr_energy_2h2p_for_1h1p_b(k_part,i_hole) += contrib + + ! Two particle correlation energy + corr_energy_2h2p_bb_2_orb(i_hole,j_hole) += contrib + corr_energy_2h2p_bb_2_orb(j_hole,i_hole) += contrib + corr_energy_2h2p_bb_2_orb(i_hole,k_part) += contrib + corr_energy_2h2p_bb_2_orb(k_part,i_hole) += contrib + corr_energy_2h2p_bb_2_orb(k_part,l_part) += contrib + corr_energy_2h2p_bb_2_orb(l_part,k_part) += contrib + + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, corr_energy_2h1p_ab_bb_per_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_a, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_b, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_double, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_2h1p] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_2h1p = 0.d0 + corr_energy_2h1p_per_orb_ab = 0.d0 + corr_energy_2h1p_per_orb_aa = 0.d0 + corr_energy_2h1p_per_orb_bb = 0.d0 + corr_energy_2h1p_ab_bb_per_2_orb = 0.d0 + corr_energy_2h1p_for_1h1p_a = 0.d0 + corr_energy_2h1p_for_1h1p_b = 0.d0 + corr_energy_2h1p_for_1h1p_double = 0.d0 + do i = 1, n_inact_orb + i_hole = list_inact(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h1p += contrib + corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib + corr_energy_2h1p_per_orb_ab(i_hole) += contrib + corr_energy_2h1p_per_orb_ab(l_part) += contrib + enddo + enddo + enddo + enddo + + ! Alpha Alpha spin correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = 1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_2h1p += contrib + corr_energy_2h1p_per_orb_aa(i_hole) += contrib + corr_energy_2h1p_per_orb_aa(l_part) += contrib + enddo + enddo + enddo + enddo + + ! Beta Beta correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = 1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib + + total_corr_e_2h1p += contrib + corr_energy_2h1p_per_orb_bb(i_hole) += contrib + corr_energy_2h1p_per_orb_aa(l_part) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_two_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_1h2p] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_1h2p = 0.d0 + corr_energy_1h2p_per_orb_ab = 0.d0 + corr_energy_1h2p_per_orb_aa = 0.d0 + corr_energy_1h2p_per_orb_bb = 0.d0 + do i = 1, n_virt_orb + i_hole = list_virt(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_ab(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib + enddo + enddo + enddo + enddo + + ! Alpha Alpha correlation energy + do i = 1, n_virt_orb + i_hole = list_virt(i) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = i+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_aa(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib + enddo + enddo + enddo + enddo + + ! Beta Beta correlation energy + do i = 1, n_virt_orb + i_hole = list_virt(i) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = i+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_bb(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, corr_energy_1h1p_spin_flip_per_orb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_1h1p_spin_flip] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_1h1p_spin_flip = 0.d0 + corr_energy_1h1p_spin_flip_per_orb = 0.d0 + do i = 1, n_inact_orb + i_hole = list_inact(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_act_orb + j_hole = list_act(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_1h1p_spin_flip += contrib + corr_energy_1h1p_spin_flip_per_orb(i_hole) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f index a4c6b652..83955e61 100644 --- a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f +++ b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f @@ -3,6 +3,7 @@ subroutine diag_inactive_virt_and_update_mos integer :: i,j,i_inact,j_inact,i_virt,j_virt double precision :: tmp(mo_tot_num_align,mo_tot_num) character*(64) :: label + print*,'Diagonalizing the occ and virt Fock operator' tmp = 0.d0 do i = 1, mo_tot_num tmp(i,i) = Fock_matrix_mo(i,i) @@ -33,3 +34,50 @@ subroutine diag_inactive_virt_and_update_mos end + +subroutine diag_inactive_virt_new_and_update_mos + implicit none + integer :: i,j,i_inact,j_inact,i_virt,j_virt,k,k_act + double precision :: tmp(mo_tot_num_align,mo_tot_num),accu,get_mo_bielec_integral_schwartz + character*(64) :: label + tmp = 0.d0 + do i = 1, mo_tot_num + tmp(i,i) = Fock_matrix_mo(i,i) + enddo + + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = i+1, n_inact_orb + j_inact = list_inact(j) + accu =0.d0 + do k = 1, n_act_orb + k_act = list_act(k) + accu += get_mo_bielec_integral_schwartz(i_inact,k_act,j_inact,k_act,mo_integrals_map) + accu -= get_mo_bielec_integral_schwartz(i_inact,k_act,k_act,j_inact,mo_integrals_map) + enddo + tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact) + accu + tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact) + accu + enddo + enddo + + do i = 1, n_virt_orb + i_virt = list_virt(i) + do j = i+1, n_virt_orb + j_virt = list_virt(j) + accu =0.d0 + do k = 1, n_act_orb + k_act = list_act(k) + accu += get_mo_bielec_integral_schwartz(i_virt,k_act,j_virt,k_act,mo_integrals_map) + enddo + tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt) - accu + tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt) - accu + enddo + enddo + + + label = "Canonical" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + soft_touch mo_coef + + +end diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index 2f662f4d..99566a8e 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -58,24 +58,24 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen call i_h_j(det_buffer(1,1,i),det_buffer(1,1,i),Nint,haa) f = 1.d0/(E_ref-haa) - if(second_order_h)then +! if(second_order_h)then lambda_i = f - else - ! You write the new Hamiltonian matrix - do k = 1, Ndet_generators - H_matrix_tmp(k,Ndet_generators+1) = H_array(k) - H_matrix_tmp(Ndet_generators+1,k) = H_array(k) - enddo - H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa - ! Then diagonalize it - call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1) - ! Then you extract the effective denominator - accu = 0.d0 - do k = 1, Ndet_generators - accu += eigenvectors(k,1) * H_array(k) - enddo - lambda_i = eigenvectors(Ndet_generators+1,1)/accu - endif +! else +! ! You write the new Hamiltonian matrix +! do k = 1, Ndet_generators +! H_matrix_tmp(k,Ndet_generators+1) = H_array(k) +! H_matrix_tmp(Ndet_generators+1,k) = H_array(k) +! enddo +! H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa +! ! Then diagonalize it +! call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1) +! ! Then you extract the effective denominator +! accu = 0.d0 +! do k = 1, Ndet_generators +! accu += eigenvectors(k,1) * H_array(k) +! enddo +! lambda_i = eigenvectors(Ndet_generators+1,1)/accu +! endif do k=1,idx(0) contrib = H_array(idx(k)) * H_array(idx(k)) * lambda_i delta_ij_generators_(idx(k), idx(k)) += contrib @@ -85,33 +85,6 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen delta_ij_generators_(idx(j), idx(k)) += contrib enddo enddo -! H_matrix_tmp_bis(idx(k),idx(k)) += contrib -! H_matrix_tmp_bis(idx(k),idx(j)) += contrib -! H_matrix_tmp_bis(idx(j),idx(k)) += contrib -! do k = 1, Ndet_generators -! do j = 1, Ndet_generators -! H_matrix_tmp_bis(k,j) = H_matrix_tmp(k,j) -! enddo -! enddo -! double precision :: H_matrix_tmp_bis(Ndet_generators,Ndet_generators) -! double precision :: eigenvectors_bis(Ndet_generators,Ndet_generators), eigenvalues_bis(Ndet_generators) -! call lapack_diag(eigenvalues_bis,eigenvectors_bis,H_matrix_tmp_bis,Ndet_generators,Ndet_generators) -! print*,'f,lambda_i = ',f,lambda_i -! print*,'eigenvalues_bi(1)',eigenvalues_bis(1) -! print*,'eigenvalues ',eigenvalues(1) -! do k = 1, Ndet_generators -! print*,'coef,coef_dres = ', eigenvectors(k,1), eigenvectors_bis(k,1) -! enddo -! pause -! accu = 0.d0 -! do k = 1, Ndet_generators -! do j = 1, Ndet_generators -! accu += eigenvectors(k,1) * eigenvectors(j,1) * (H_matrix_tmp(k,j) + delta_ij_generators_(k,j)) -! enddo -! enddo -! print*,'accu,eigv = ',accu,eigenvalues(1) -! pause - enddo end diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f new file mode 100644 index 00000000..8656b633 --- /dev/null +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -0,0 +1,59 @@ +program foboscf + implicit none + call run_prepare + no_oa_or_av_opt = .True. + touch no_oa_or_av_opt + call routine_fobo_scf + call save_mos + +end + +subroutine run_prepare + implicit none + no_oa_or_av_opt = .False. + touch no_oa_or_av_opt + call damping_SCF + call diag_inactive_virt_and_update_mos +end + +subroutine routine_fobo_scf + implicit none + integer :: i,j + print*,'' + print*,'' + character*(64) :: label + label = "Natural" + do i = 1, 5 + print*,'*******************************************************************************' + print*,'*******************************************************************************' + print*,'FOBO-SCF Iteration ',i + print*,'*******************************************************************************' + print*,'*******************************************************************************' + if(speed_up_convergence_foboscf)then + if(i==3)then + threshold_lmct = max(threshold_lmct,0.001) + threshold_mlct = max(threshold_mlct,0.05) + soft_touch threshold_lmct threshold_mlct + endif + if(i==4)then + threshold_lmct = max(threshold_lmct,0.005) + threshold_mlct = max(threshold_mlct,0.07) + soft_touch threshold_lmct threshold_mlct + endif + if(i==5)then + threshold_lmct = max(threshold_lmct,0.01) + threshold_mlct = max(threshold_mlct,0.1) + soft_touch threshold_lmct threshold_mlct + endif + endif + call FOBOCI_lmct_mlct_old_thr + call save_osoci_natural_mos + call damping_SCF + call diag_inactive_virt_and_update_mos + call clear_mo_map + call provide_properties + enddo + + + +end diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index 087f791b..dc6519b8 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -9,12 +9,9 @@ subroutine FOBOCI_lmct_mlct_old_thr double precision :: norm_tmp(N_states),norm_total(N_states) logical :: test_sym double precision :: thr,hij - double precision :: threshold double precision, allocatable :: dressing_matrix(:,:) logical :: verbose,is_ok verbose = .True. - threshold = threshold_singles - print*,'threshold = ',threshold thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -36,7 +33,14 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'' print*,'' print*,'DOING FIRST LMCT !!' + print*,'Threshold_lmct = ',threshold_lmct + integer(bit_kind) , allocatable :: zero_bitmask(:,:) + integer(bit_kind) , allocatable :: psi_singles(:,:,:) + logical :: lmct + double precision, allocatable :: psi_singles_coef(:,:) + allocate( zero_bitmask(N_int,2) ) do i = 1, n_inact_orb + lmct = .True. integer :: i_hole_osoci i_hole_osoci = list_inact(i) print*,'--------------------------' @@ -51,27 +55,91 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'Passed set generators' call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) - call is_a_good_candidate(threshold,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle - ! so all the mono excitation on the new generators allocate(dressing_matrix(N_det_generators,N_det_generators)) + dressing_matrix = 0.d0 if(.not.do_it_perturbative)then -! call all_single - dressing_matrix = 0.d0 + do k = 1, N_det_generators do l = 1, N_det_generators call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) dressing_matrix(k,l) = hkl enddo enddo - double precision :: hkl -! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) -! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) - call debug_det(reunion_of_bitmask,N_int) + hkl = dressing_matrix(1,1) + do k = 1, N_det_generators + dressing_matrix(k,k) = dressing_matrix(k,k) - hkl + enddo + print*,'Naked matrix' + do k = 1, N_det_generators + write(*,'(100(F12.5,X))')dressing_matrix(k,:) + enddo + + ! Do all the single excitations on top of the CAS and 1h determinants + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single +! if(dressing_2h2p)then +! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct) +! endif + +! ! Change the mask of the holes and particles to perform all the +! ! double excitations that starts from the active space in order +! ! to introduce the Coulomb hole in the active space +! ! These are the 1h2p excitations that have the i_hole_osoci hole in common +! ! and the 2p if there is more than one electron in the active space +! do k = 1, N_int +! zero_bitmask(k,1) = 0_bit_kind +! zero_bitmask(k,2) = 0_bit_kind +! enddo +! ! hole is possible only in the orbital i_hole_osoci +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int) +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int) +! ! and in the active space +! do k = 1, n_act_orb +! call set_bit_to_integer(list_act(k),zero_bitmask(1,1),N_int) +! call set_bit_to_integer(list_act(k),zero_bitmask(1,2),N_int) +! enddo +! call set_bitmask_hole_as_input(zero_bitmask) + +! call set_bitmask_particl_as_input(reunion_of_bitmask) + +! call all_1h2p +! call diagonalize_CI_SC2 +! call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) + +! ! Change the mask of the holes and particles to perform all the +! ! double excitations that from the orbital i_hole_osoci +! do k = 1, N_int +! zero_bitmask(k,1) = 0_bit_kind +! zero_bitmask(k,2) = 0_bit_kind +! enddo +! ! hole is possible only in the orbital i_hole_osoci +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int) +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int) +! call set_bitmask_hole_as_input(zero_bitmask) + +! call set_bitmask_particl_as_input(reunion_of_bitmask) + +! call set_psi_det_to_generators +! call all_2h2p +! call diagonalize_CI_SC2 + double precision :: hkl + call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) + hkl = dressing_matrix(1,1) + do k = 1, N_det_generators + dressing_matrix(k,k) = dressing_matrix(k,k) - hkl + enddo + print*,'Dressed matrix' + do k = 1, N_det_generators + write(*,'(100(F12.5,X))')dressing_matrix(k,:) + enddo +! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) endif call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) + do k = 1, N_states print*,'norm_tmp = ',norm_tmp(k) norm_total(k) += norm_tmp(k) @@ -83,9 +151,12 @@ subroutine FOBOCI_lmct_mlct_old_thr if(.True.)then print*,'' print*,'DOING THEN THE MLCT !!' + print*,'Threshold_mlct = ',threshold_mlct + lmct = .False. do i = 1, n_virt_orb integer :: i_particl_osoci i_particl_osoci = list_virt(i) + print*,'--------------------------' ! First set the current generators to the one of restart call set_generators_to_generators_restart @@ -107,7 +178,7 @@ subroutine FOBOCI_lmct_mlct_old_thr call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) !! ! so all the mono excitation on the new generators - call is_a_good_candidate(threshold,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle allocate(dressing_matrix(N_det_generators,N_det_generators)) @@ -122,6 +193,9 @@ subroutine FOBOCI_lmct_mlct_old_thr ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) call all_single +! if(dressing_2h2p)then +! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct) +! endif endif call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) do k = 1, N_states @@ -132,24 +206,6 @@ subroutine FOBOCI_lmct_mlct_old_thr deallocate(dressing_matrix) enddo endif - if(.False.)then - print*,'LAST loop for all the 1h-1p' - print*,'--------------------------' - ! First set the current generators to the one of restart - call set_generators_to_generators_restart - call set_psi_det_to_generators - call initialize_bitmask_to_restart_ones - ! Impose that only the hole i_hole_osoci can be done - call set_bitmask_particl_as_input(inact_virt_bitmask) - call set_bitmask_hole_as_input(inact_virt_bitmask) -! call set_bitmask_particl_as_input(reunion_of_bitmask) -! call set_bitmask_hole_as_input(reunion_of_bitmask) - call all_single - call set_intermediate_normalization_1h1p(norm_tmp) - norm_total += norm_tmp - call update_density_matrix_osoci - endif - print*,'norm_total = ',norm_total norm_total = norm_generators_restart @@ -174,10 +230,8 @@ subroutine FOBOCI_mlct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - double precision :: threshold logical :: verbose,is_ok verbose = .False. - threshold = 1.d-2 thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -216,7 +270,7 @@ subroutine FOBOCI_mlct_old call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) ! ! so all the mono excitation on the new generators - call is_a_good_candidate(threshold,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose) print*,'is_ok = ',is_ok is_ok =.True. if(.not.is_ok)cycle @@ -250,10 +304,8 @@ subroutine FOBOCI_lmct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - double precision :: threshold logical :: verbose,is_ok verbose = .False. - threshold = 1.d-2 thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -290,7 +342,7 @@ subroutine FOBOCI_lmct_old call set_generators_to_psi_det call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) - call is_a_good_candidate(threshold,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle ! ! so all the mono excitation on the new generators diff --git a/plugins/FOBOCI/foboci_reunion.irp.f b/plugins/FOBOCI/foboci_reunion.irp.f new file mode 100644 index 00000000..fcfaff58 --- /dev/null +++ b/plugins/FOBOCI/foboci_reunion.irp.f @@ -0,0 +1,18 @@ +program osoci_program +implicit none + do_it_perturbative = .True. + touch do_it_perturbative + call FOBOCI_lmct_mlct_old_thr + call provide_all_the_rest +end +subroutine provide_all_the_rest +implicit none +integer :: i + call update_one_body_dm_mo + call set_lmct_mlct_to_psi_det + call diagonalize_CI + call save_wavefunction + + +end + diff --git a/plugins/FOBOCI/generators_restart_save.irp.f b/plugins/FOBOCI/generators_restart_save.irp.f index dca4c901..09d4aa2b 100644 --- a/plugins/FOBOCI/generators_restart_save.irp.f +++ b/plugins/FOBOCI/generators_restart_save.irp.f @@ -1,126 +1,74 @@ -use bitmasks +use bitmasks + BEGIN_PROVIDER [ integer, N_det_generators_restart ] implicit none BEGIN_DOC - ! Number of determinants in the wave function + ! Read the wave function END_DOC - logical :: exists - character*64 :: label + integer :: i integer, save :: ifirst = 0 -!if(ifirst == 0)then - PROVIDE ezfio_filename - call ezfio_has_determinants_n_det(exists) - print*,'exists = ',exists - if(.not.exists)then - print*,'The OSOCI needs a restart WF' - print*,'There are none in the EZFIO file ...' - print*,'Stopping ...' - stop - endif - print*,'passed N_det_generators_restart' - call ezfio_get_determinants_n_det(N_det_generators_restart) - ASSERT (N_det_generators_restart > 0) + double precision :: norm + if(ifirst == 0)then + call ezfio_get_determinants_n_det(N_det_generators_restart) ifirst = 1 -!endif + else + print*,'PB in generators_restart restart !!!' + endif + call write_int(output_determinants,N_det_generators_restart,'Number of generators_restart') END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,N_det_generators_restart) ] &BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ] implicit none BEGIN_DOC - ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file - ! is empty + ! read wf + ! END_DOC - integer :: i - logical :: exists - character*64 :: label - + integer :: i, k integer, save :: ifirst = 0 -!if(ifirst == 0)then - provide N_det_generators_restart - if(.True.)then - call ezfio_has_determinants_N_int(exists) - if (exists) then - call ezfio_has_determinants_bit_kind(exists) - if (exists) then - call ezfio_has_determinants_N_det(exists) - if (exists) then - call ezfio_has_determinants_N_states(exists) - if (exists) then - call ezfio_has_determinants_psi_det(exists) - endif - endif - endif - endif - - if(.not.exists)then - print*,'The OSOCI needs a restart WF' - print*,'There are none in the EZFIO file ...' - print*,'Stopping ...' - stop - endif - print*,'passed psi_det_generators_restart' - - call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) - do i = 1, N_int - ref_generators_restart(i,1) = psi_det_generators_restart(i,1,1) - ref_generators_restart(i,2) = psi_det_generators_restart(i,2,1) - enddo - endif + double precision, allocatable :: psi_coef_read(:,:) + if(ifirst == 0)then + call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) + do k = 1, N_int + ref_generators_restart(k,1) = psi_det_generators_restart(k,1,1) + ref_generators_restart(k,2) = psi_det_generators_restart(k,2,1) + enddo + allocate (psi_coef_read(N_det_generators_restart,N_states)) + call ezfio_get_determinants_psi_coef(psi_coef_read) + do k = 1, N_states + do i = 1, N_det_generators_restart + psi_coef_generators_restart(i,k) = psi_coef_read(i,k) + enddo + enddo ifirst = 1 -!endif + deallocate(psi_coef_read) + else + print*,'PB in generators_restart restart !!!' + endif END_PROVIDER - - -BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (psi_det_size,N_states_diag) ] - implicit none - BEGIN_DOC - ! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file - ! is empty - END_DOC - - integer :: i,k, N_int2 - logical :: exists - double precision, allocatable :: psi_coef_read(:,:) - character*(64) :: label - - integer, save :: ifirst = 0 -!if(ifirst == 0)then - psi_coef_generators_restart = 0.d0 - do i=1,N_states_diag - psi_coef_generators_restart(i,i) = 1.d0 - enddo - - call ezfio_has_determinants_psi_coef(exists) - - if(.not.exists)then - print*,'The OSOCI needs a restart WF' - print*,'There are none in the EZFIO file ...' - print*,'Stopping ...' - stop - endif - print*,'passed psi_coef_generators_restart' - - if (exists) then - - allocate (psi_coef_read(N_det_generators_restart,N_states)) - call ezfio_get_determinants_psi_coef(psi_coef_read) - do k=1,N_states - do i=1,N_det_generators_restart - psi_coef_generators_restart(i,k) = psi_coef_read(i,k) - enddo - enddo - deallocate(psi_coef_read) - - endif - ifirst = 1 -!endif - - - +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 END_PROVIDER +BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] + implicit none + BEGIN_DOC + ! Memo to skip useless selectors + END_DOC + select_max = huge(1.d0) +END_PROVIDER + + BEGIN_PROVIDER [ integer, N_det_generators ] +&BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,10000) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (10000,N_states) ] + +END_PROVIDER diff --git a/plugins/FOBOCI/hcc_1h1p.irp.f b/plugins/FOBOCI/hcc_1h1p.irp.f new file mode 100644 index 00000000..66cf2fd4 --- /dev/null +++ b/plugins/FOBOCI/hcc_1h1p.irp.f @@ -0,0 +1,83 @@ +program test_sc2 + implicit none + read_wf = .True. + touch read_wf + call routine + + +end + +subroutine routine + implicit none + double precision, allocatable :: energies(:),diag_H_elements(:) + double precision, allocatable :: H_matrix(:,:) + allocate(energies(N_states),diag_H_elements(N_det)) + call diagonalize_CI + call test_hcc + call test_mulliken +! call SC2_1h1p(psi_det,psi_coef,energies, & +! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + allocate(H_matrix(N_det,N_det)) + call SC2_1h1p_full(psi_det,psi_coef,energies, & + H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + deallocate(H_matrix) + integer :: i,j + double precision :: accu,coef_hf +! coef_hf = 1.d0/psi_coef(1,1) +! do i = 1, N_det +! psi_coef(i,1) *= coef_hf +! enddo + touch psi_coef + call pouet +end + +subroutine pouet + implicit none + double precision :: accu,coef_hf +! provide one_body_dm_mo_alpha one_body_dm_mo_beta +! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int) +! touch one_body_dm_mo_alpha one_body_dm_mo_beta + call test_hcc + call test_mulliken +! call save_wavefunction + +end + +subroutine test_hcc + implicit none + double precision :: accu + integer :: i,j + print*,'Z AU GAUSS MHZ cm^-1' + do i = 1, nucl_num + write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) + enddo + +end + +subroutine test_mulliken + double precision :: accu + integer :: i + integer :: j + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),mulliken_spin_densities(i) + accu += mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu +!print*,'AO SPIN POPULATIONS' + accu = 0.d0 +!do i = 1, ao_num +! accu += spin_gross_orbital_product(i) +! write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) +!enddo +!print*,'sum = ',accu +!accu = 0.d0 +!print*,'Angular momentum analysis' +!do i = 0, ao_l_max +! accu += spin_population_angular_momentum(i) +! print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) +!print*,'sum = ',accu +!enddo + +end + diff --git a/plugins/FOBOCI/modify_generators.irp.f b/plugins/FOBOCI/modify_generators.irp.f index c756f0c2..359b6405 100644 --- a/plugins/FOBOCI/modify_generators.irp.f +++ b/plugins/FOBOCI/modify_generators.irp.f @@ -6,6 +6,7 @@ subroutine set_generators_to_psi_det END_DOC N_det_generators = N_det integer :: i,k + print*,'N_det = ',N_det do i=1,N_det_generators do k=1,N_int psi_det_generators(k,1,i) = psi_det(k,1,i) diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 49dcafc3..8e2f2e53 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -24,6 +24,7 @@ subroutine new_approach double precision, allocatable :: dressing_matrix_1h1p(:,:) double precision, allocatable :: dressing_matrix_2h1p(:,:) double precision, allocatable :: dressing_matrix_1h2p(:,:) + double precision, allocatable :: dressing_matrix_extra_1h_or_1p(:,:) double precision, allocatable :: H_matrix_tmp(:,:) logical :: verbose,is_ok @@ -45,7 +46,7 @@ subroutine new_approach verbose = .True. - threshold = threshold_singles + threshold = threshold_lmct print*,'threshold = ',threshold thr = 1.d-12 print*,'' @@ -81,12 +82,14 @@ subroutine new_approach ! so all the mono excitation on the new generators allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) allocate(dressing_matrix_2h1p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)) dressing_matrix_1h1p = 0.d0 dressing_matrix_2h1p = 0.d0 + dressing_matrix_extra_1h_or_1p = 0.d0 if(.not.do_it_perturbative)then n_good_hole +=1 ! call all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) - call all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + call all_single_for_1h(i_hole_foboci,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p) allocate(H_matrix_tmp(N_det_generators,N_det_generators)) do j = 1,N_det_generators do k = 1, N_det_generators @@ -96,7 +99,7 @@ subroutine new_approach enddo do j = 1, N_det_generators do k = 1, N_det_generators - H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + dressing_matrix_extra_1h_or_1p(j,k) enddo enddo hjk = H_matrix_tmp(1,1) @@ -130,6 +133,7 @@ subroutine new_approach endif deallocate(dressing_matrix_1h1p) deallocate(dressing_matrix_2h1p) + deallocate(dressing_matrix_extra_1h_or_1p) enddo print*,'' @@ -155,12 +159,14 @@ subroutine new_approach ! so all the mono excitation on the new generators allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) allocate(dressing_matrix_1h2p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)) dressing_matrix_1h1p = 0.d0 dressing_matrix_1h2p = 0.d0 + dressing_matrix_extra_1h_or_1p = 0.d0 if(.not.do_it_perturbative)then n_good_hole +=1 ! call all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) - call all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + call all_single_for_1p(i_particl_osoci,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) allocate(H_matrix_tmp(N_det_generators,N_det_generators)) do j = 1,N_det_generators do k = 1, N_det_generators @@ -170,7 +176,7 @@ subroutine new_approach enddo do j = 1, N_det_generators do k = 1, N_det_generators - H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + dressing_matrix_extra_1h_or_1p(j,k) enddo enddo hjk = H_matrix_tmp(1,1) @@ -205,7 +211,10 @@ subroutine new_approach endif deallocate(dressing_matrix_1h1p) deallocate(dressing_matrix_1h2p) + deallocate(dressing_matrix_extra_1h_or_1p) enddo + + double precision, allocatable :: H_matrix_total(:,:) integer :: n_det_total n_det_total = N_det_generators_restart + n_good_det @@ -221,7 +230,7 @@ subroutine new_approach !!! Adding the averaged dressing coming from the 1h1p that are redundant for each of the "n_good_hole" 1h H_matrix_total(i,j) += dressing_matrix_restart_1h1p(i,j)/dble(n_good_hole+n_good_particl) !!! Adding the dressing coming from the 2h1p that are not redundant for the any of CI calculations - H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + dressing_matrix_restart_1h2p(i,j) enddo enddo do i = 1, n_good_det @@ -244,25 +253,79 @@ subroutine new_approach H_matrix_total(n_det_generators_restart+j,n_det_generators_restart+i) = hij enddo enddo - print*,'H matrix to diagonalize' - double precision :: href - href = H_matrix_total(1,1) - do i = 1, n_det_total - H_matrix_total(i,i) -= href + + ! Adding the correlation energy + logical :: orb_taken_good_det(mo_tot_num) + double precision :: phase + integer :: n_h,n_p,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + integer :: h1,h2,p1,p2,s1,s2 + logical, allocatable :: one_hole_or_one_p(:) + integer, allocatable :: holes_or_particle(:) + allocate(one_hole_or_one_p(n_good_det), holes_or_particle(n_good_det)) + orb_taken_good_det = .False. + do i = 1, n_good_det + n_h = number_of_holes(psi_good_det(1,1,i)) + n_p = number_of_particles(psi_good_det(1,1,i)) + call get_excitation(ref_bitmask,psi_good_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(n_h == 0 .and. n_p == 1)then + orb_taken_good_det(h1) = .True. + one_hole_or_one_p(i) = .True. + holes_or_particle(i) = h1 + endif + if(n_h == 1 .and. n_p == 0)then + orb_taken_good_det(p1) = .True. + one_hole_or_one_p(i) = .False. + holes_or_particle(i) = p1 + endif enddo - do i = 1, n_det_total - write(*,'(100(X,F16.8))')H_matrix_total(i,:) - enddo - double precision, allocatable :: eigvalues(:),eigvectors(:,:) - allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total)) - call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total) - print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion + href - do i = 1, n_det_total - print*,'coef = ',eigvectors(i,1) - enddo - integer(bit_kind), allocatable :: psi_det_final(:,:,:) - double precision, allocatable :: psi_coef_final(:,:) - double precision :: norm + + do i = 1, N_det_generators_restart + ! Add the 2h2p, 2h1p and 1h2p correlation energy + H_matrix_total(i,i) += total_corr_e_2h2p + total_corr_e_2h1p + total_corr_e_1h2p + total_corr_e_1h1p_spin_flip + ! Substract the 2h1p part that have already been taken into account + do j = 1, n_inact_orb + iorb = list_inact(j) + if(.not.orb_taken_good_det(iorb))cycle + H_matrix_total(i,i) -= corr_energy_2h1p_per_orb_ab(iorb) - corr_energy_2h1p_per_orb_bb(iorb) - corr_energy_1h1p_spin_flip_per_orb(iorb) + enddo + ! Substract the 1h2p part that have already been taken into account + do j = 1, n_virt_orb + iorb = list_virt(j) + if(.not.orb_taken_good_det(iorb))cycle + H_matrix_total(i,i) -= corr_energy_1h2p_per_orb_ab(iorb) - corr_energy_1h2p_per_orb_aa(iorb) + enddo + enddo + + do i = 1, N_good_det + ! Repeat the 2h2p correlation energy + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_2h2p + ! Substract the part that can not be repeated + ! If it is a 1h + if(one_hole_or_one_p(i))then + ! 2h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h2p_per_orb_bb(holes_or_particle(i)) + ! You can repeat a certain part of the 1h2p correlation energy + ! that is everything except the part that involves the hole of the 1h + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_1h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_1h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_1h2p_per_orb_bb(holes_or_particle(i)) + + else + ! 2h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h2p_per_orb_aa(holes_or_particle(i)) + ! You can repeat a certain part of the 2h1p correlation energy + ! that is everything except the part that involves the hole of the 1p + ! 2h1p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h1p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h1p_per_orb_aa(holes_or_particle(i)) + endif + enddo + allocate(psi_coef_final(n_det_total, N_states)) allocate(psi_det_final(N_int,2,n_det_total)) do i = 1, N_det_generators_restart @@ -277,22 +340,222 @@ subroutine new_approach psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i) enddo enddo - norm = 0.d0 + + + double precision :: href + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + integer(bit_kind), allocatable :: psi_det_final(:,:,:) + double precision, allocatable :: psi_coef_final(:,:) + double precision :: norm + allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total)) + + call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total) + print*,'' + print*,'' + print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1) + print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion do i = 1, n_det_total - do j = 1, N_states - psi_coef_final(i,j) = eigvectors(i,j) - enddo - norm += psi_coef_final(i,1)**2 -! call debug_det(psi_det_final(1, 1, i), N_int) + print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1) enddo - print*,'norm = ',norm + + integer(bit_kind), allocatable :: psi_det_remaining_1h_or_1p(:,:,:) + integer(bit_kind), allocatable :: key_tmp(:,:) + integer :: n_det_remaining_1h_or_1p + integer :: ispin,i_ok + allocate(key_tmp(N_int,2),psi_det_remaining_1h_or_1p(N_int,2,n_inact_orb*n_act_orb+n_virt_orb*n_act_orb)) + logical :: is_already_present + logical, allocatable :: one_hole_or_one_p_bis(:) + integer, allocatable :: holes_or_particle_bis(:) + double precision,allocatable :: H_array(:) + allocate(one_hole_or_one_p_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb), holes_or_particle_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb)) + allocate(H_array(n_det_total)) + ! Dressing with the remaining 1h determinants + print*,'' + print*,'' + print*,'Dressing with the remaining 1h determinants' + n_det_remaining_1h_or_1p = 0 + do i = 1, n_inact_orb + iorb = list_inact(i) + if(orb_taken_good_det(iorb))cycle + do j = 1, n_act_orb + jorb = list_act(j) + ispin = 2 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,iorb,jorb,ispin,i_ok) + if(i_ok .ne.1)cycle + is_already_present = .False. + H_array = 0.d0 + call i_h_j(key_tmp,key_tmp,N_int,hij) + href = ref_bitmask_energy - hij + href = 1.d0/href + do k = 1, n_det_total + call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int) + if(degree == 0)then + is_already_present = .True. + exit + endif + enddo + if(is_already_present)cycle + n_det_remaining_1h_or_1p +=1 + one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .True. + holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb + do k = 1, N_int + psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1) + psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2) + enddo + ! do k = 1, n_det_total + ! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij) + ! H_array(k) = hij + ! enddo + ! do k = 1, n_det_total + ! do l = 1, n_det_total + ! H_matrix_total(k,l) += H_array(k) * H_array(l) * href + ! enddo + ! enddo + enddo + enddo + ! Dressing with the remaining 1p determinants + print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p + print*,'Dressing with the remaining 1p determinants' + do i = 1, n_virt_orb + iorb = list_virt(i) + if(orb_taken_good_det(iorb))cycle + do j = 1, n_act_orb + jorb = list_act(j) + ispin = 1 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,jorb,iorb,ispin,i_ok) + if(i_ok .ne.1)cycle + is_already_present = .False. + H_array = 0.d0 + call i_h_j(key_tmp,key_tmp,N_int,hij) + href = ref_bitmask_energy - hij + href = 1.d0/href + do k = 1, n_det_total + call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int) + if(degree == 0)then + is_already_present = .True. + exit + endif + enddo + if(is_already_present)cycle + n_det_remaining_1h_or_1p +=1 + one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .False. + holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb + do k = 1, N_int + psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1) + psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2) + enddo +! do k = 1, n_det_total +! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij) +! H_array(k) = hij +! enddo +! do k = 1, n_det_total +! do l = 1, n_det_total +! H_matrix_total(k,l) += H_array(k) * H_array(l) * href +! enddo +! enddo + enddo + enddo + print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p + deallocate(key_tmp,H_array) + + double precision, allocatable :: eigvalues_bis(:),eigvectors_bis(:,:),H_matrix_total_bis(:,:) + integer :: n_det_final + n_det_final = n_det_total + n_det_remaining_1h_or_1p + allocate(eigvalues_bis(n_det_final),eigvectors_bis(n_det_final,n_det_final),H_matrix_total_bis(n_det_final,n_det_final)) + print*,'passed the allocate, building the big matrix' + do i = 1, n_det_total + do j = 1, n_det_total + H_matrix_total_bis(i,j) = H_matrix_total(i,j) + enddo + enddo + do i = 1, n_det_remaining_1h_or_1p + do j = 1, n_det_remaining_1h_or_1p + call i_h_j(psi_det_remaining_1h_or_1p(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij) + H_matrix_total_bis(n_det_total+i,n_det_total+j) = hij + enddo + enddo + do i = 1, n_det_total + do j = 1, n_det_remaining_1h_or_1p + call i_h_j(psi_det_final(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij) + H_matrix_total_bis(i,n_det_total+j) = hij + H_matrix_total_bis(n_det_total+j,i) = hij + enddo + enddo + print*,'passed the matrix' + do i = 1, n_det_remaining_1h_or_1p + if(one_hole_or_one_p_bis(i))then + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h2p_per_orb_bb(holes_or_particle_bis(i)) + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_1h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_1h2p_per_orb_bb(holes_or_particle_bis(i)) + else + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h2p_per_orb_aa(holes_or_particle_bis(i)) + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_2h1p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h1p_per_orb_aa(holes_or_particle_bis(i)) + + endif + enddo + do i = 2, n_det_final + do j = i+1, n_det_final + H_matrix_total_bis(i,j) = 0.d0 + H_matrix_total_bis(j,i) = 0.d0 + enddo + enddo + do i = 1, n_det_final + write(*,'(500(F10.5,X))')H_matrix_total_bis(i,:) + enddo + call lapack_diag(eigvalues_bis,eigvectors_bis,H_matrix_total_bis,n_det_final,n_det_final) + print*,'e_dressed = ',eigvalues_bis(1) + nuclear_repulsion + do i = 1, n_det_final + print*,'coef = ',eigvectors_bis(i,1),H_matrix_total_bis(i,i) - H_matrix_total_bis(1,1) + enddo + do j = 1, N_states + do i = 1, n_det_total + psi_coef_final(i,j) = eigvectors_bis(i,j) + norm += psi_coef_final(i,j)**2 + enddo + norm = 1.d0/dsqrt(norm) + do i = 1, n_det_total + psi_coef_final(i,j) = psi_coef_final(i,j) * norm + enddo + enddo + + + deallocate(eigvalues_bis,eigvectors_bis,H_matrix_total_bis) + + +!print*,'H matrix to diagonalize' +!href = H_matrix_total(1,1) +!do i = 1, n_det_total +! H_matrix_total(i,i) -= href +!enddo +!do i = 1, n_det_total +! write(*,'(100(X,F16.8))')H_matrix_total(i,:) +!enddo +!call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total) +!print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1) +!print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion +!do i = 1, n_det_total +! print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1) +!enddo +!norm = 0.d0 +!do i = 1, n_det_total +! do j = 1, N_states +! psi_coef_final(i,j) = eigvectors(i,j) +! enddo +! norm += psi_coef_final(i,1)**2 +!enddo +!print*,'norm = ',norm call set_psi_det_as_input_psi(n_det_total,psi_det_final,psi_coef_final) - print*,'' -!do i = 1, N_det -! call debug_det(psi_det(1,1,i),N_int) -! print*,'coef = ',psi_coef(i,1) -!enddo + + do i = 1, N_det + call debug_det(psi_det(1,1,i),N_int) + print*,'coef = ',psi_coef(i,1) + enddo provide one_body_dm_mo integer :: i_core,iorb,jorb,i_inact,j_inact,i_virt,j_virt,j_core @@ -360,14 +623,14 @@ subroutine new_approach print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) endif diff --git a/plugins/FOBOCI/new_new_approach.irp.f b/plugins/FOBOCI/new_new_approach.irp.f new file mode 100644 index 00000000..b904a5b3 --- /dev/null +++ b/plugins/FOBOCI/new_new_approach.irp.f @@ -0,0 +1,132 @@ +program test_new_new + implicit none + read_wf = .True. + touch read_wf + call test +end + + +subroutine test + implicit none + integer :: i,j,k,l + call diagonalize_CI + call set_generators_to_psi_det + print*,'Initial coefficients' + do i = 1, N_det + print*,'' + call debug_det(psi_det(1,1,i),N_int) + print*,'psi_coef = ',psi_coef(i,1) + print*,'' + enddo + double precision, allocatable :: dressing_matrix(:,:) + double precision :: hij + double precision :: phase + integer :: n_h,n_p,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + integer :: h1,h2,p1,p2,s1,s2 + allocate(dressing_matrix(N_det_generators,N_det_generators)) + do i = 1, N_det_generators + do j = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij) + dressing_matrix(i,j) = hij + enddo + enddo + href = dressing_matrix(1,1) + print*,'Diagonal part of the dressing' + do i = 1, N_det_generators + print*,'delta e = ',dressing_matrix(i,i) - href + enddo + call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) + double precision :: href + print*,'' + ! One considers that the following excitation classes are not repeatable on the 1h and 1p determinants : + ! + 1h1p spin flip + ! + 2h1p + ! + 1h2p + ! But the 2h2p are correctly taken into account +!dressing_matrix(1,1) += total_corr_e_1h2p + total_corr_e_2h1p + total_corr_e_1h1p_spin_flip +!do i = 1, N_det_generators +! dressing_matrix(i,i) += total_corr_e_2h2p +! n_h = number_of_holes(psi_det(1,1,i)) +! n_p = number_of_particles(psi_det(1,1,i)) +! if(n_h == 1 .and. n_p ==0)then +! +! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! print*,'' +! print*,' 1h det ' +! print*,'' +! call debug_det(psi_det_generators(1,1,i),N_int) +! print*,'h1,p1 = ',h1,p1 +! print*,'total_corr_e_2h2p ',total_corr_e_2h2p +! print*,'corr_energy_2h2p_per_orb_ab(h1)',corr_energy_2h2p_per_orb_ab(h1) +! print*,'corr_energy_2h2p_per_orb_bb(h1)',corr_energy_2h2p_per_orb_bb(h1) +! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1) +! dressing_matrix(1,1) += -corr_energy_2h1p_per_orb_aa(h1) - corr_energy_2h1p_per_orb_ab(h1) -corr_energy_2h1p_per_orb_bb(h1) & +! -corr_energy_1h1p_spin_flip_per_orb(h1) +! endif +! if(n_h == 0 .and. n_p ==1)then +! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! print*,'' +! print*,' 1p det ' +! print*,'' +! call debug_det(psi_det_generators(1,1,i),N_int) +! print*,'h1,p1 = ',h1,p1 +! print*,'total_corr_e_2h2p ',total_corr_e_2h2p +! print*,'corr_energy_2h2p_per_orb_ab(p1)',corr_energy_2h2p_per_orb_ab(p1) +! print*,'corr_energy_2h2p_per_orb_aa(p1)',corr_energy_2h2p_per_orb_aa(p1) +! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(p1) - corr_energy_2h2p_per_orb_aa(p1) +! dressing_matrix(1,1) += -corr_energy_1h2p_per_orb_aa(p1) - corr_energy_1h2p_per_orb_ab(p1) -corr_energy_1h2p_per_orb_bb(p1) +! endif +!enddo +!href = dressing_matrix(1,1) +!print*,'Diagonal part of the dressing' +!do i = 1, N_det_generators +! print*,'delta e = ',dressing_matrix(i,i) - href +!enddo + call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) + print*,'After dressing matrix' + print*,'' + print*,'' + do i = 1, N_det + print*,'psi_coef = ',psi_coef(i,1) + enddo +!print*,'' +!print*,'' +!print*,'Canceling the dressing part of the interaction between 1h and 1p' +!do i = 2, N_det_generators +! do j = i+1, N_det_generators +! call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij) +! dressing_matrix(i,j) = hij +! dressing_matrix(j,i) = hij +! enddo +!enddo +!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) +!print*,'' +!print*,'' +!do i = 1, N_det +! print*,'psi_coef = ',psi_coef(i,1) +!enddo +!print*,'' +!print*,'' +!print*,'Canceling the interaction between 1h and 1p' + +!print*,'' +!print*,'' +!do i = 2, N_det_generators +! do j = i+1, N_det_generators +! dressing_matrix(i,j) = 0.d0 +! dressing_matrix(j,i) = 0.d0 +! enddo +!enddo +!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) +!do i = 1, N_det +! print*,'psi_coef = ',psi_coef(i,1) +!enddo + call save_natural_mos + deallocate(dressing_matrix) + + +end diff --git a/plugins/FOBOCI/routines_dressing.irp.f b/plugins/FOBOCI/routines_dressing.irp.f index 910f1109..125143da 100644 --- a/plugins/FOBOCI/routines_dressing.irp.f +++ b/plugins/FOBOCI/routines_dressing.irp.f @@ -55,15 +55,11 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det i_pert = 0 endif do j = 1, ndet_generators_input - if(dabs(H_array(j)*lambda_i).gt.0.5d0)then + if(dabs(H_array(j)*lambda_i).gt.0.1d0)then i_pert = 1 exit endif enddo -! print*,'' -! print*,'lambda_i,f = ',lambda_i,f -! print*,'i_pert = ',i_pert -! print*,'' if(i_pert==1)then lambda_i = f i_pert_count +=1 @@ -79,9 +75,122 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det enddo enddo enddo + href = dressing_matrix(1,1) + print*,'Diagonal part of the dressing' + do i = 1, ndet_generators_input + print*,'delta e = ',dressing_matrix(i,i) - href + enddo !print*,'i_pert_count = ',i_pert_count end +subroutine update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,H_jj_in) + use bitmasks + implicit none + integer, intent(in) :: ndet_generators_input + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input) + double precision, intent(in) :: H_jj_in(N_det) + double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input) + integer :: i,j,n_det_ref_tmp,degree + double precision :: href + n_det_ref_tmp = 0 + do i = 1, N_det + do j = 1, Ndet_generators_input + call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int) + if(degree == 0)then + dressing_matrix(j,j) += H_jj_in(i) + n_det_ref_tmp +=1 + exit + endif + enddo + enddo + if( ndet_generators_input .ne. n_det_ref_tmp)then + print*,'Problem !!!! ' + print*,' ndet_generators .ne. n_det_ref_tmp !!!' + print*,'ndet_generators,n_det_ref_tmp' + print*,ndet_generators_input,n_det_ref_tmp + stop + endif + + href = dressing_matrix(1,1) + print*,'' + print*,'Update with the SC2 dressing' + print*,'' + print*,'Diagonal part of the dressing' + do i = 1, ndet_generators_input + print*,'delta e = ',dressing_matrix(i,i) - href + enddo +end + +subroutine provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, & + psi_det_outer_input,psi_coef_outer_input,n_det_outer_input) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_input + integer(bit_kind), intent(in) :: psi_det_ref_input(N_int,2,n_det_ref_input) + double precision, intent(in) :: psi_coef_ref_input(n_det_ref_input,N_states) + integer, intent(in) :: n_det_outer_input + integer(bit_kind), intent(in) :: psi_det_outer_input(N_int,2,n_det_outer_input) + double precision, intent(in) :: psi_coef_outer_input(n_det_outer_input,N_states) + + double precision, intent(inout) :: dressing_matrix(n_det_ref_input,n_det_ref_input) + + + integer :: i_pert, i_pert_count,i,j,k + double precision :: f,href,hka,lambda_i + double precision :: H_array(n_det_ref_input),accu + integer :: n_h_out,n_p_out,n_p_in,n_h_in,number_of_holes,number_of_particles + call i_h_j(psi_det_ref_input(1,1,1),psi_det_ref_input(1,1,1),N_int,href) + i_pert_count = 0 + do i = 1, n_det_outer_input + call i_h_j(psi_det_outer_input(1,1,i),psi_det_outer_input(1,1,i),N_int,hka) + f = 1.d0/(href - hka) + H_array = 0.d0 + accu = 0.d0 +! n_h_out = number_of_holes(psi_det_outer_input(1,1,i)) +! n_p_out = number_of_particles(psi_det_outer_input(1,1,i)) + do j=1,n_det_ref_input + n_h_in = number_of_holes(psi_det_ref_input(1,1,j)) + n_p_in = number_of_particles(psi_det_ref_input(1,1,j)) +! if(n_h_in == 0 .and. n_h_in == 0)then + call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka) +! else +! hka = 0.d0 +! endif + H_array(j) = hka + accu += psi_coef_ref_input(j,1) * hka + enddo + lambda_i = psi_coef_outer_input(i,1)/accu + i_pert = 1 + if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then + i_pert = 0 + endif + do j = 1, n_det_ref_input + if(dabs(H_array(j)*lambda_i).gt.0.5d0)then + i_pert = 1 + exit + endif + enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! i_pert = 0 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if(i_pert==1)then + lambda_i = f + i_pert_count +=1 + endif + do k=1,n_det_ref_input + double precision :: contrib + contrib = H_array(k) * H_array(k) * lambda_i + dressing_matrix(k, k) += contrib + do j=k+1,n_det_ref_input + contrib = H_array(k) * H_array(j) * lambda_i + dressing_matrix(k, j) += contrib + dressing_matrix(j, k) += contrib + enddo + enddo + enddo +end + + subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, & psi_det_outer_input,psi_coef_outer_input,n_det_outer_input) use bitmasks @@ -112,16 +221,17 @@ subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi accu += psi_coef_ref_input(j,1) * hka enddo lambda_i = psi_coef_outer_input(i,1)/accu - i_pert = 1 + i_pert = 0 if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then i_pert = 0 endif do j = 1, n_det_ref_input - if(dabs(H_array(j)*lambda_i).gt.0.3d0)then + if(dabs(H_array(j)*lambda_i).gt.0.5d0)then i_pert = 1 exit endif enddo +! i_pert = 0 if(i_pert==1)then lambda_i = f i_pert_count +=1 @@ -170,114 +280,379 @@ subroutine diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_ end -subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) +subroutine give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p) use bitmasks implicit none - integer, intent(out) :: n_det_1h1p, n_det_2h1p + integer, intent(in) :: i_hole + integer, intent(out) :: n_det_1h1p, n_det_2h1p,n_det_extra_1h_or_1p integer :: i integer :: n_det_ref_restart_tmp,n_det_1h integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det n_det_ref_restart_tmp = 0 n_det_1h = 0 n_det_1h1p = 0 n_det_2h1p = 0 + n_det_extra_1h_or_1p = 0 do i = 1, N_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) if(n_h == 0 .and. n_p == 0)then n_det_ref_restart_tmp +=1 else if (n_h ==1 .and. n_p==0)then - n_det_1h +=1 + if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then + n_det_1h +=1 + else + n_det_extra_1h_or_1p +=1 + endif + else if (n_h ==0 .and. n_p==1)then + n_det_extra_1h_or_1p +=1 else if (n_h ==1 .and. n_p==1)then n_det_1h1p +=1 else if (n_h ==2 .and. n_p==1)then n_det_2h1p +=1 else print*,'PB !!!!' - print*,'You have something else than a 1h, 1h1p or 2h1p' + print*,'You have something else than a 1h, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p call debug_det(psi_det(1,1,i),N_int) stop endif enddo -! if(n_det_1h.ne.1)then -! print*,'PB !! You have more than one 1h' -! stop -! endif if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then print*,'PB !!!!' print*,'You have forgotten something in your generators ... ' stop endif - + if(n_det_2h1p + n_det_1h1p + n_det_extra_1h_or_1p + n_det_generators .ne. N_det)then + print*,'PB !!!!' + print*,'You have forgotten something in your generators ... ' + stop + endif end -subroutine give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) +subroutine give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p) use bitmasks implicit none - integer, intent(out) :: n_det_1h1p, n_det_1h2p + integer, intent(out) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p integer :: i integer :: n_det_ref_restart_tmp,n_det_1h integer :: number_of_holes,n_h, number_of_particles,n_p - n_det_ref_restart_tmp = 0 - n_det_1h = 0 + logical :: is_the_hole_in_det + n_det_ref_1h_1p = 0 + n_det_2h1p = 0 n_det_1h1p = 0 - n_det_1h2p = 0 do i = 1, N_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) if(n_h == 0 .and. n_p == 0)then - n_det_ref_restart_tmp +=1 + n_det_ref_1h_1p +=1 + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p +=1 else if (n_h ==0 .and. n_p==1)then - n_det_1h +=1 + n_det_ref_1h_1p +=1 + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p +=1 + else if (n_h ==2 .and. n_p==1)then + n_det_2h1p +=1 + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + +end + +subroutine give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p) + use bitmasks + implicit none + integer, intent(out) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + n_det_ref_1h_1p = 0 + n_det_1h2p = 0 + n_det_1h1p = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_1h_1p +=1 + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p +=1 + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p +=1 else if (n_h ==1 .and. n_p==1)then n_det_1h1p +=1 else if (n_h ==1 .and. n_p==2)then n_det_1h2p +=1 else print*,'PB !!!!' - print*,'You have something else than a 1p, 1h1p or 1h2p' + print*,'You have something else than a 1h, 1p, 1h1p or 1h2p' + print*,'n_h,n_p = ',n_h,n_p call debug_det(psi_det(1,1,i),N_int) stop endif enddo - if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then + +end + +subroutine give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p + integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p) + integer(bit_kind), intent(out) :: psi_det_2h1p(N_int,2,n_det_2h1p) + integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p) + double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states) + double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states) + integer :: n_det_ref_1h_1p_tmp,n_det_2h1p_tmp,n_det_1h1p_tmp + integer :: i,j + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + integer, allocatable :: index_ref_1h_1p(:) + integer, allocatable :: index_2h1p(:) + integer, allocatable :: index_1h1p(:) + allocate(index_ref_1h_1p(n_det)) + allocate(index_2h1p(n_det)) + allocate(index_1h1p(n_det)) + n_det_ref_1h_1p_tmp = 0 + n_det_2h1p_tmp = 0 + n_det_1h1p_tmp = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p_tmp +=1 + index_1h1p(n_det_1h1p_tmp) = i + else if (n_h ==2 .and. n_p==1)then + n_det_2h1p_tmp +=1 + index_2h1p(n_det_2h1p_tmp) = i + else print*,'PB !!!!' - print*,'You have forgotten something in your generators ... ' - stop - endif + print*,'You have something else than a 1h, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + do i = 1, n_det_2h1p + do j = 1, N_int + psi_det_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i)) + psi_det_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i)) + enddo + do j = 1, N_states + psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j) + enddo + enddo + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i)) + enddo + do j = 1, N_states + psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j) + enddo + enddo + + do i = 1, n_det_ref_1h_1p + do j = 1, N_int + psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i)) + psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i)) + enddo + do j = 1, N_states + psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j) + enddo + enddo + +end + +subroutine give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p) + integer(bit_kind), intent(out) :: psi_det_1h2p(N_int,2,n_det_1h2p) + integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p) + double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states) + double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states) + integer :: n_det_ref_1h_1p_tmp,n_det_1h2p_tmp,n_det_1h1p_tmp + integer :: i,j + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + integer, allocatable :: index_ref_1h_1p(:) + integer, allocatable :: index_1h2p(:) + integer, allocatable :: index_1h1p(:) + allocate(index_ref_1h_1p(n_det)) + allocate(index_1h2p(n_det)) + allocate(index_1h1p(n_det)) + n_det_ref_1h_1p_tmp = 0 + n_det_1h2p_tmp = 0 + n_det_1h1p_tmp = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p_tmp +=1 + index_1h1p(n_det_1h1p_tmp) = i + else if (n_h ==1 .and. n_p==2)then + n_det_1h2p_tmp +=1 + index_1h2p(n_det_1h2p_tmp) = i + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1p, 1h1p or 1h2p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + do i = 1, n_det_1h2p + do j = 1, N_int + psi_det_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i)) + psi_det_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i)) + enddo + do j = 1, N_states + psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j) + enddo + enddo + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i)) + enddo + do j = 1, N_states + psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j) + enddo + enddo + + do i = 1, n_det_ref_1h_1p + do j = 1, N_int + psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i)) + psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i)) + enddo + do j = 1, N_states + psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j) + enddo + enddo + +end + + + +subroutine give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p) + use bitmasks + implicit none + integer, intent(in) ::i_particl + integer, intent(out) :: n_det_1h1p, n_det_1h2p,n_det_extra_1h_or_1p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1p + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_particl_in_det + n_det_ref_restart_tmp = 0 + n_det_1p = 0 + n_det_1h1p = 0 + n_det_1h2p = 0 + n_det_extra_1h_or_1p = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_restart_tmp +=1 + else if (n_h ==0 .and. n_p==1)then + if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then + n_det_1p +=1 + else + n_det_extra_1h_or_1p +=1 + endif + else if (n_h ==1 .and. n_p==0)then + n_det_extra_1h_or_1p +=1 + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p +=1 + else if (n_h ==1 .and. n_p==2)then + n_det_1h2p +=1 + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1p, 1h1p or 1h2p' + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo +!if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then +! print*,'PB !!!!' +! print*,'You have forgotten something in your generators ... ' +! stop +!endif end -subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) +subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) use bitmasks implicit none - integer, intent(in) :: n_det_1h1p,n_det_2h1p + integer, intent(in) :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p,i_hole integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) integer(bit_kind), intent(out) :: psi_2h1p(N_int,2,n_det_2h1p) + integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p) double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p, N_states) + double precision, intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states) integer :: i,j integer :: degree integer :: number_of_holes,n_h, number_of_particles,n_p - integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp,n_det_extra_1h_or_1p_tmp + integer :: n_det_1h_tmp integer, allocatable :: index_generator(:) integer, allocatable :: index_1h1p(:) integer, allocatable :: index_2h1p(:) + integer, allocatable :: index_extra_1h_or_1p(:) + logical :: is_the_hole_in_det allocate(index_1h1p(n_det)) allocate(index_2h1p(n_det)) + allocate(index_extra_1h_or_1p(n_det)) allocate(index_generator(N_det)) n_det_generators_tmp = 0 n_det_1h1p_tmp = 0 n_det_2h1p_tmp = 0 + n_det_extra_1h_or_1p_tmp = 0 + n_det_1h_tmp = 0 do i = 1, n_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) @@ -287,6 +662,16 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o else if (n_h ==2 .and. n_p==1)then n_det_2h1p_tmp +=1 index_2h1p(n_det_2h1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then + n_det_1h_tmp +=1 + else + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + endif endif do j = 1, N_det_generators call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) @@ -315,6 +700,12 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o stop endif + if(n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp)then + print*,'PB !!!' + print*,'n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp' + stop + endif + do i = 1,N_det_generators do j = 1, N_int psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i)) @@ -345,41 +736,59 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o enddo enddo + do i = 1, n_det_extra_1h_or_1p + do j = 1, N_int + psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i)) + psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i)) + enddo + do j = 1, N_states + psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j) + enddo + enddo deallocate(index_generator) deallocate(index_1h1p) deallocate(index_2h1p) + deallocate(index_extra_1h_or_1p) end -subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) +subroutine split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) use bitmasks implicit none - integer, intent(in) :: n_det_1h1p,n_det_1h2p + integer, intent(in) :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p,i_particl integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) integer(bit_kind), intent(out) :: psi_1h2p(N_int,2,n_det_1h2p) + integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p) double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p, N_states) + double precision, intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states) integer :: i,j integer :: degree integer :: number_of_holes,n_h, number_of_particles,n_p - integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp,n_det_extra_1h_or_1p_tmp integer, allocatable :: index_generator(:) integer, allocatable :: index_1h1p(:) integer, allocatable :: index_1h2p(:) + integer, allocatable :: index_extra_1h_or_1p(:) + logical :: is_the_particl_in_det + integer :: n_det_1p_tmp allocate(index_1h1p(n_det)) allocate(index_1h2p(n_det)) + allocate(index_extra_1h_or_1p(n_det)) allocate(index_generator(N_det)) n_det_generators_tmp = 0 n_det_1h1p_tmp = 0 n_det_1h2p_tmp = 0 + n_det_extra_1h_or_1p_tmp = 0 + n_det_1p_tmp = 0 do i = 1, n_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) @@ -389,6 +798,15 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o else if (n_h ==1 .and. n_p==2)then n_det_1h2p_tmp +=1 index_1h2p(n_det_1h2p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then + n_det_1p_tmp +=1 + else + n_det_extra_1h_or_1p_tmp +=1 + endif endif do j = 1, N_det_generators call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) @@ -448,9 +866,20 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o enddo + do i = 1, n_det_extra_1h_or_1p + do j = 1, N_int + psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i)) + psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i)) + enddo + do j = 1, N_states + psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j) + enddo + enddo + deallocate(index_generator) deallocate(index_1h1p) deallocate(index_1h2p) + deallocate(index_extra_1h_or_1p) end diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 696011a9..4aca60d7 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -332,20 +332,20 @@ subroutine save_osoci_natural_mos enddo tmp = tmp_bis -!! Symetrization act-virt - do j = 1, n_virt_orb - j_virt= list_virt(j) - accu = 0.d0 - do i = 1, n_act_orb - jorb = list_act(i) - accu += dabs(tmp_bis(j_virt,jorb)) - enddo - do i = 1, n_act_orb - iorb = list_act(i) - tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) - tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) - enddo - enddo +!!! Symetrization act-virt +! do j = 1, n_virt_orb +! j_virt= list_virt(j) +! accu = 0.d0 +! do i = 1, n_act_orb +! jorb = list_act(i) +! accu += dabs(tmp_bis(j_virt,jorb)) +! enddo +! do i = 1, n_act_orb +! iorb = list_act(i) +! tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) +! tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) +! enddo +! enddo !! Symetrization act-inact !do j = 1, n_inact_orb @@ -387,16 +387,16 @@ subroutine save_osoci_natural_mos print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo enddo @@ -410,8 +410,9 @@ subroutine save_osoci_natural_mos enddo label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) - soft_touch mo_coef +!soft_touch mo_coef deallocate(tmp,occ) @@ -518,16 +519,16 @@ subroutine set_osoci_natural_mos print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo enddo @@ -602,15 +603,210 @@ end subroutine provide_properties implicit none - integer :: i - double precision :: accu - if(.True.)then - accu= 0.d0 - do i = 1, nucl_num - accu += mulliken_spin_densities(i) - print*,i,nucl_charge(i),mulliken_spin_densities(i) - enddo - print*,'Sum of Mulliken SD = ',accu - endif + call print_mulliken_sd + call print_hcc end + + + subroutine dress_diag_elem_2h1p(dressing_H_mat_elem,ndet,lmct,i_hole) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet,i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 0.and.n_p==0)then ! CAS + dressing_H_mat_elem(i)+= total_corr_e_2h1p + if(lmct)then + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(i_hole) - corr_energy_2h1p_per_orb_bb(i_hole) + endif + endif + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(p1) - corr_energy_2h1p_per_orb_aa(p1) + else if (n_h == 1.and.n_p==1)then ! 1h1p +! if(degree==1)then + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i)+= - corr_energy_2h1p_per_orb_ab(h1) +! else +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) +! dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) +! endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(h1) + endif + enddo + + end + + subroutine dress_diag_elem_1h2p(dressing_H_mat_elem,ndet,lmct,i_hole) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet,i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 0.and.n_p==0)then ! CAS + dressing_H_mat_elem(i)+= total_corr_e_1h2p + if(.not.lmct)then + dressing_H_mat_elem(i) += - corr_energy_1h2p_per_orb_ab(i_hole) - corr_energy_1h2p_per_orb_aa(i_hole) + endif + endif + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i)+= total_corr_e_1h2p - corr_energy_1h2p_per_orb_ab(h1) + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 1.and.n_p==1)then ! 1h1p + if(degree==1)then + dressing_H_mat_elem(i)+= total_corr_e_1h2p + dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) + else + dressing_H_mat_elem(i) +=0.d0 + endif +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) +! dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) +! endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i)+= total_corr_e_1h2p + dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) - corr_energy_1h2p_per_orb_ab(h1) + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i) += 0.d0 + endif + enddo + + end + + subroutine dress_diag_elem_2h2p(dressing_H_mat_elem,ndet) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + dressing_H_mat_elem(i)+= total_corr_e_2h2p + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + else if (n_h == 1.and.n_p==1)then ! 1h1p + if(degree==1)then + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_a(h1,p1) + corr_energy_2h2p_for_1h1p_b(h1,p1)) + else + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) + dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) + endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1) & + - corr_energy_2h2p_per_orb_ab(h2) & + - 0.5d0 * ( corr_energy_2h2p_per_orb_bb(h2) + corr_energy_2h2p_per_orb_bb(h2)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) + if(s1.ne.s2)then + dressing_H_mat_elem(i) += corr_energy_2h2p_ab_2_orb(h1,h2) + else + dressing_H_mat_elem(i) += corr_energy_2h2p_bb_2_orb(h1,h2) + endif + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) + if(s1.ne.s2)then + dressing_H_mat_elem(i) += corr_energy_2h2p_ab_2_orb(p1,p2) + else + dressing_H_mat_elem(i) += corr_energy_2h2p_bb_2_orb(p1,p2) + endif + endif + enddo + + end + + subroutine diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole,lmct) + implicit none + double precision, allocatable :: dressing_H_mat_elem(:),energies(:) + integer, intent(in) :: i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + integer :: i + double precision :: hij + allocate(dressing_H_mat_elem(N_det),energies(N_states_diag)) + print*,'' + print*,'dressing with the 2h2p in a CC logic' + print*,'' + do i = 1, N_det + call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij) + dressing_H_mat_elem(i) = hij + enddo + call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det) + call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,i_hole) + call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,i_hole) + call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states_diag,N_int,output_determinants) + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo + + + deallocate(dressing_H_mat_elem) + + + + end diff --git a/plugins/Generators_restart/generators.irp.f b/plugins/Generators_restart/generators.irp.f index 0a82e6f9..17854330 100644 --- a/plugins/Generators_restart/generators.irp.f +++ b/plugins/Generators_restart/generators.irp.f @@ -1,5 +1,5 @@ use bitmasks - + BEGIN_PROVIDER [ integer, N_det_generators ] implicit none BEGIN_DOC @@ -8,17 +8,18 @@ BEGIN_PROVIDER [ integer, N_det_generators ] integer :: i integer, save :: ifirst = 0 double precision :: norm - read_wf = .True. if(ifirst == 0)then - N_det_generators = N_det + call ezfio_get_determinants_n_det(N_det_generators) ifirst = 1 + else + print*,'PB in generators restart !!!' endif call write_int(output_determinants,N_det_generators,'Number of generators') END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ] implicit none BEGIN_DOC ! read wf @@ -26,17 +27,20 @@ END_PROVIDER END_DOC integer :: i, k integer, save :: ifirst = 0 + double precision, allocatable :: psi_coef_read(:,:) if(ifirst == 0)then - do i=1,N_det_generators - do k=1,N_int - psi_det_generators(k,1,i) = psi_det(k,1,i) - psi_det_generators(k,2,i) = psi_det(k,2,i) - enddo + call read_dets(psi_det_generators,N_int,N_det_generators) + allocate (psi_coef_read(N_det_generators,N_states)) + call ezfio_get_determinants_psi_coef(psi_coef_read) do k = 1, N_states - psi_coef_generators(i,k) = psi_coef(i,k) + do i = 1, N_det_generators + psi_coef_generators(i,k) = psi_coef_read(i,k) + enddo enddo - enddo ifirst = 1 + deallocate(psi_coef_read) + else + print*,'PB in generators restart !!!' endif END_PROVIDER diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index a407860f..aa6f02b0 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,9 @@ subroutine damping_SCF write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') '====','================','================','================', '====' write(output_hartree_fock,*) - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + if(.not.no_oa_or_av_opt)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + endif call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') call ezfio_set_hartree_fock_energy(E_min) diff --git a/plugins/Molden/NEEDED_CHILDREN_MODULES b/plugins/Molden/NEEDED_CHILDREN_MODULES index 305dfb78..80d0af12 100644 --- a/plugins/Molden/NEEDED_CHILDREN_MODULES +++ b/plugins/Molden/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -MO_Basis Utils +MO_Basis Utils diff --git a/plugins/Molden/aos.irp.f b/plugins/Molden/aos.irp.f deleted file mode 100644 index 71f8c5b8..00000000 --- a/plugins/Molden/aos.irp.f +++ /dev/null @@ -1,196 +0,0 @@ -BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ] - implicit none - BEGIN_DOC -! ao_l = l value of the AO: a+b+c in x^a y^b z^c - END_DOC - integer :: i - do i=1,ao_num - ao_l_char(i) = l_to_character(ao_l(i)) - enddo -END_PROVIDER - - -BEGIN_PROVIDER [ character*(128), l_to_character, (0:4)] - BEGIN_DOC - ! character corresponding to the "L" value of an AO orbital - END_DOC - implicit none - l_to_character(0)='S' - l_to_character(1)='P' - l_to_character(2)='D' - l_to_character(3)='F' - l_to_character(4)='G' -END_PROVIDER - - BEGIN_PROVIDER [ integer, Nucl_N_Aos, (nucl_num)] -&BEGIN_PROVIDER [ integer, N_AOs_max ] - implicit none - integer :: i - BEGIN_DOC - ! Number of AOs per atom - END_DOC - Nucl_N_Aos = 0 - do i = 1, ao_num - Nucl_N_Aos(ao_nucl(i)) +=1 - enddo - N_AOs_max = maxval(Nucl_N_Aos) -END_PROVIDER - - BEGIN_PROVIDER [ integer, Nucl_Aos, (nucl_num,N_AOs_max)] - implicit none - BEGIN_DOC - ! List of AOs attached on each atom - END_DOC - integer :: i - integer, allocatable :: nucl_tmp(:) - allocate(nucl_tmp(nucl_num)) - nucl_tmp = 0 - Nucl_Aos = 0 - do i = 1, ao_num - nucl_tmp(ao_nucl(i))+=1 - Nucl_Aos(ao_nucl(i),nucl_tmp(ao_nucl(i))) = i - enddo - deallocate(nucl_tmp) -END_PROVIDER - - - BEGIN_PROVIDER [ integer, Nucl_list_shell_Aos, (nucl_num,N_AOs_max)] -&BEGIN_PROVIDER [ integer, Nucl_num_shell_Aos, (nucl_num)] - implicit none - integer :: i,j,k - BEGIN_DOC - ! Index of the shell type Aos and of the corresponding Aos - ! Per convention, for P,D,F and G AOs, we take the index - ! of the AO with the the corresponding power in the "X" axis - END_DOC - do i = 1, nucl_num - Nucl_num_shell_Aos(i) = 0 - - do j = 1, Nucl_N_Aos(i) - if(ao_l(Nucl_Aos(i,j))==0)then - ! S type function - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - elseif(ao_l(Nucl_Aos(i,j))==1)then - ! P type function - if(ao_power(Nucl_Aos(i,j),1)==1)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - elseif(ao_l(Nucl_Aos(i,j))==2)then - ! D type function - if(ao_power(Nucl_Aos(i,j),1)==2)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - elseif(ao_l(Nucl_Aos(i,j))==3)then - ! F type function - if(ao_power(Nucl_Aos(i,j),1)==3)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - elseif(ao_l(Nucl_Aos(i,j))==4)then - ! G type function - if(ao_power(Nucl_Aos(i,j),1)==4)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - endif - - enddo - enddo - -END_PROVIDER - - -BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] - implicit none - integer :: i - character*(4) :: give_ao_character_space - do i=1,ao_num - - if(ao_l(i)==0)then - ! S type AO - give_ao_character_space = 'S ' - elseif(ao_l(i) == 1)then - ! P type AO - if(ao_power(i,1)==1)then - give_ao_character_space = 'X ' - elseif(ao_power(i,2) == 1)then - give_ao_character_space = 'Y ' - else - give_ao_character_space = 'Z ' - endif - elseif(ao_l(i) == 2)then - ! D type AO - if(ao_power(i,1)==2)then - give_ao_character_space = 'XX ' - elseif(ao_power(i,2) == 2)then - give_ao_character_space = 'YY ' - elseif(ao_power(i,3) == 2)then - give_ao_character_space = 'ZZ ' - elseif(ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'XY ' - elseif(ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'XZ ' - else - give_ao_character_space = 'YZ ' - endif - elseif(ao_l(i) == 3)then - ! F type AO - if(ao_power(i,1)==3)then - give_ao_character_space = 'XXX ' - elseif(ao_power(i,2) == 3)then - give_ao_character_space = 'YYY ' - elseif(ao_power(i,3) == 3)then - give_ao_character_space = 'ZZZ ' - elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'XXY ' - elseif(ao_power(i,1) == 2 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'XXZ ' - elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1)then - give_ao_character_space = 'YYX ' - elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'YYZ ' - elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1)then - give_ao_character_space = 'ZZX ' - elseif(ao_power(i,3) == 2 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'ZZY ' - elseif(ao_power(i,3) == 1 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'XYZ ' - endif - elseif(ao_l(i) == 4)then - ! G type AO - if(ao_power(i,1)==4)then - give_ao_character_space = 'XXXX' - elseif(ao_power(i,2) == 4)then - give_ao_character_space = 'YYYY' - elseif(ao_power(i,3) == 4)then - give_ao_character_space = 'ZZZZ' - elseif(ao_power(i,1) == 3 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'XXXY' - elseif(ao_power(i,1) == 3 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'XXXZ' - elseif(ao_power(i,2) == 3 .and. ao_power(i,1) == 1)then - give_ao_character_space = 'YYYX' - elseif(ao_power(i,2) == 3 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'YYYZ' - elseif(ao_power(i,3) == 3 .and. ao_power(i,1) == 1)then - give_ao_character_space = 'ZZZX' - elseif(ao_power(i,3) == 3 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'ZZZY' - elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 2)then - give_ao_character_space = 'XXYY' - elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 2)then - give_ao_character_space = 'YYZZ' - elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'XXYZ' - elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then - give_ao_character_space = 'YYXZ' - elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then - give_ao_character_space = 'ZZXY' - endif - endif - ao_l_char_space(i) = give_ao_character_space - enddo -END_PROVIDER diff --git a/plugins/Molden/print_mo.irp.f b/plugins/Molden/print_mo.irp.f index b147fe50..6ac51bdb 100644 --- a/plugins/Molden/print_mo.irp.f +++ b/plugins/Molden/print_mo.irp.f @@ -104,6 +104,8 @@ subroutine write_Ao_basis(i_unit_output) write(i_unit_output,*)'' write(i_unit_output,'(A47,2X,I3)')'TOTAL NUMBER OF BASIS SET SHELLS =', i_shell write(i_unit_output,'(A47,2X,I3)')'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =', ao_num +! this is for the new version of molden + write(i_unit_output,'(A12)')'PP =NONE' write(i_unit_output,*)'' @@ -126,7 +128,9 @@ subroutine write_Mo_basis(i_unit_output) write(i_unit_output,'(18X,F8.5)')-1.d0 write(i_unit_output,*)'' do i = 1, ao_num - write(i_unit_output,'(2X,I3, 2X A1, I3, 2X A4 , F9.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j) +! write(i_unit_output,'(2X,I3, 2X A1, I3, 2X A4 , F9.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j) +! F12.6 for larger coefficients... + write(i_unit_output,'(2X,I3, 2X A1, I3, 2X A4 , F12.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j) ! write(i_unit_output,'(I3, X A1, X I3, X A4 X F16.8)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)) enddo write(i_unit_output,*)'' diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index e990a37c..e406cd03 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -125,6 +125,8 @@ subroutine pt2_moller_plesset ($arguments) delta_e = (Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)) + & (Fock_matrix_diag_mo(h2) - Fock_matrix_diag_mo(p2)) delta_e = 1.d0/delta_e +! print*,'h1,p1',h1,p1 +! print*,'h2,p2',h2,p2 else if (degree == 1) then call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1) diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f index c1d88d2c..e31b3ba4 100644 --- a/plugins/Properties/hyperfine_constants.irp.f +++ b/plugins/Properties/hyperfine_constants.irp.f @@ -133,3 +133,16 @@ END_PROVIDER enddo END_PROVIDER + + +subroutine print_hcc + implicit none + double precision :: accu + integer :: i,j + print*,'Z AU GAUSS MHZ cm^-1' + do i = 1, nucl_num + write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) + enddo + +end + diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f index d56c9a44..cc0a2f8e 100644 --- a/plugins/Properties/mulliken.irp.f +++ b/plugins/Properties/mulliken.irp.f @@ -105,3 +105,34 @@ END_PROVIDER enddo END_PROVIDER + + +subroutine print_mulliken_sd + implicit none + double precision :: accu + integer :: i + integer :: j + print*,'Mulliken spin densities' + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),mulliken_spin_densities(i) + accu += mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu + print*,'AO SPIN POPULATIONS' + accu = 0.d0 + do i = 1, ao_num + accu += spin_gross_orbital_product(i) + write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) + enddo + print*,'sum = ',accu + accu = 0.d0 + print*,'Angular momentum analysis' + do i = 0, ao_l_max + accu += spin_population_angular_momentum(i) + print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) + print*,'sum = ',accu + enddo + +end + diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f index f0091e1e..45bca5e6 100644 --- a/plugins/Properties/print_hcc.irp.f +++ b/plugins/Properties/print_hcc.irp.f @@ -1,17 +1,6 @@ -program print_hcc +program print_hcc_main implicit none read_wf = .True. touch read_wf - call test + call print_hcc end -subroutine test - implicit none - double precision :: accu - integer :: i,j - print*,'Z AU GAUSS MHZ cm^-1' - do i = 1, nucl_num - write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) - enddo - -end - diff --git a/plugins/Properties/print_mulliken.irp.f b/plugins/Properties/print_mulliken.irp.f index 100c8556..d4be534a 100644 --- a/plugins/Properties/print_mulliken.irp.f +++ b/plugins/Properties/print_mulliken.irp.f @@ -2,34 +2,5 @@ program print_mulliken implicit none read_wf = .True. touch read_wf - print*,'Mulliken spin densities' - - call test + call print_mulliken_sd end -subroutine test - double precision :: accu - integer :: i - integer :: j - accu= 0.d0 - do i = 1, nucl_num - print*,i,nucl_charge(i),mulliken_spin_densities(i) - accu += mulliken_spin_densities(i) - enddo - print*,'Sum of Mulliken SD = ',accu - print*,'AO SPIN POPULATIONS' - accu = 0.d0 - do i = 1, ao_num - accu += spin_gross_orbital_product(i) - write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) - enddo - print*,'sum = ',accu - accu = 0.d0 - print*,'Angular momentum analysis' - do i = 0, ao_l_max - accu += spin_population_angular_momentum(i) - print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) - print*,'sum = ',accu - enddo - -end - diff --git a/plugins/loc_cele/loc.f b/plugins/loc_cele/loc.f index 575932a3..edc3aa7a 100644 --- a/plugins/loc_cele/loc.f +++ b/plugins/loc_cele/loc.f @@ -17,9 +17,11 @@ C data small/1.d-6/ zprt=.true. - niter=500 + niter=1000000 conv=1.d-8 +C niter=1000000 +C conv=1.d-6 write (6,5) n,m,conv 5 format (//5x,'Unitary transformation of',i3,' vectors'/ * 5x,'following the principle of maximum overlap with a set of', diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f index e9c26f9d..c9036aa1 100644 --- a/plugins/loc_cele/loc_cele.irp.f +++ b/plugins/loc_cele/loc_cele.irp.f @@ -92,13 +92,182 @@ - nrot(1) = 6 ! number of orbitals to be localized + nrot(1) = 64 ! number of orbitals to be localized integer :: index_rot(1000,1) cmoref = 0.d0 + irot = 0 + +! H2 molecule for the mixed localization + do i=1,64 + irot(i,1) = i+2 + enddo + + do i=1,17 + cmoref(i+1,i,1)=1.d0 + enddo + cmoref(19,19-1,1)=1.d0 + cmoref(20,19-1,1)=-1.d0 + cmoref(19,20-1,1)=-1.d0 + cmoref(20,20-1,1)=-1.d0 + cmoref(21,20-1,1)=2.d0 + cmoref(22,21-1,1)=1.d0 + cmoref(23,22-1,1)=1.d0 + cmoref(24,23-1,1)=1.d0 + + + cmoref(25,24-1,1)=1.d0 + cmoref(26,24-1,1)=-1.d0 + cmoref(25,25-1,1)=-1.d0 + cmoref(26,25-1,1)=-1.d0 + cmoref(27,25-1,1)=2.d0 + cmoref(28,26-1,1)=1.d0 + cmoref(29,27-1,1)=1.d0 + cmoref(30,28-1,1)=1.d0 + + cmoref(31,29-1,1)=1.d0 + cmoref(32,29-1,1)=-1.d0 + cmoref(31,30-1,1)=-1.d0 + cmoref(32,30-1,1)=-1.d0 + cmoref(33,30-1,1)=2.d0 + cmoref(34,31-1,1)=1.d0 + cmoref(35,32-1,1)=1.d0 + cmoref(36,33-1,1)=1.d0 + + do i=33,49 + cmoref(i+5,i,1)= 1.d0 + enddo + + cmoref(55,52-2,1)=1.d0 + cmoref(56,52-2,1)=-1.d0 + cmoref(55,53-2,1)=-1.d0 + cmoref(56,53-2,1)=-1.d0 + cmoref(57,53-2,1)=2.d0 + cmoref(58,54-2,1)=1.d0 + cmoref(59,55-2,1)=1.d0 + cmoref(60,56-2,1)=1.d0 + + cmoref(61,57-2,1)=1.d0 + cmoref(62,57-2,1)=-1.d0 + cmoref(61,58-2,1)=-1.d0 + cmoref(62,58-2,1)=-1.d0 + cmoref(63,58-2,1)=2.d0 + cmoref(64,59-2,1)=1.d0 + cmoref(65,60-2,1)=1.d0 + cmoref(66,61-2,1)=1.d0 + + cmoref(67,62-2,1)=1.d0 + cmoref(68,62-2,1)=-1.d0 + cmoref(67,63-2,1)=-1.d0 + cmoref(68,63-2,1)=-1.d0 + cmoref(69,63-2,1)=2.d0 + cmoref(70,64-2,1)=1.d0 + cmoref(71,65-2,1)=1.d0 + cmoref(72,66-2,1)=1.d0 +! H2 molecule +! do i=1,66 +! irot(i,1) = i +! enddo +! +! do i=1,18 +! cmoref(i,i,1)=1.d0 +! enddo +! cmoref(19,19,1)=1.d0 +! cmoref(20,19,1)=-1.d0 +! cmoref(19,20,1)=-1.d0 +! cmoref(20,20,1)=-1.d0 +! cmoref(21,20,1)=2.d0 +! cmoref(22,21,1)=1.d0 +! cmoref(23,22,1)=1.d0 +! cmoref(24,23,1)=1.d0 +! +! +! cmoref(25,24,1)=1.d0 +! cmoref(26,24,1)=-1.d0 +! cmoref(25,25,1)=-1.d0 +! cmoref(26,25,1)=-1.d0 +! cmoref(27,25,1)=2.d0 +! cmoref(28,26,1)=1.d0 +! cmoref(29,27,1)=1.d0 +! cmoref(30,28,1)=1.d0 +! +! cmoref(31,29,1)=1.d0 +! cmoref(32,29,1)=-1.d0 +! cmoref(31,30,1)=-1.d0 +! cmoref(32,30,1)=-1.d0 +! cmoref(33,30,1)=2.d0 +! cmoref(34,31,1)=1.d0 +! cmoref(35,32,1)=1.d0 +! cmoref(36,33,1)=1.d0 +! +! do i=34,51 +! cmoref(i+3,i,1)= 1.d0 +! enddo +! +! cmoref(55,52,1)=1.d0 +! cmoref(56,52,1)=-1.d0 +! cmoref(55,53,1)=-1.d0 +! cmoref(56,53,1)=-1.d0 +! cmoref(57,53,1)=2.d0 +! cmoref(58,54,1)=1.d0 +! cmoref(59,55,1)=1.d0 +! cmoref(60,56,1)=1.d0 +! +! cmoref(61,57,1)=1.d0 +! cmoref(62,57,1)=-1.d0 +! cmoref(61,58,1)=-1.d0 +! cmoref(62,58,1)=-1.d0 +! cmoref(63,58,1)=2.d0 +! cmoref(64,59,1)=1.d0 +! cmoref(65,60,1)=1.d0 +! cmoref(66,61,1)=1.d0 +! +! cmoref(67,62,1)=1.d0 +! cmoref(68,62,1)=-1.d0 +! cmoref(67,63,1)=-1.d0 +! cmoref(68,63,1)=-1.d0 +! cmoref(69,63,1)=2.d0 +! cmoref(70,64,1)=1.d0 +! cmoref(71,65,1)=1.d0 +! cmoref(72,66,1)=1.d0 +! H atom +! do i=1,33 +! irot(i,1) = i +! enddo +! +! do i=1,18 +! cmoref(i,i,1)=1.d0 +! enddo +! cmoref(19,19,1)=1.d0 +! cmoref(20,19,1)=-1.d0 +! cmoref(19,20,1)=-1.d0 +! cmoref(20,20,1)=-1.d0 +! cmoref(21,20,1)=2.d0 +! cmoref(22,21,1)=1.d0 +! cmoref(23,22,1)=1.d0 +! cmoref(24,23,1)=1.d0 + + +! cmoref(25,24,1)=1.d0 +! cmoref(26,24,1)=-1.d0 +! cmoref(25,25,1)=-1.d0 +! cmoref(26,25,1)=-1.d0 +! cmoref(27,25,1)=2.d0 +! cmoref(28,26,1)=1.d0 +! cmoref(29,27,1)=1.d0 +! cmoref(30,28,1)=1.d0 +! +! cmoref(31,29,1)=1.d0 +! cmoref(32,29,1)=-1.d0 +! cmoref(31,30,1)=-1.d0 +! cmoref(32,30,1)=-1.d0 +! cmoref(33,30,1)=2.d0 +! cmoref(34,31,1)=1.d0 +! cmoref(35,32,1)=1.d0 +! cmoref(36,33,1)=1.d0 ! Definition of the index of the MO to be rotated ! irot(2,1) = 21 ! the first mo to be rotated is the 21 th MO @@ -106,25 +275,67 @@ ! irot(4,1) = 23 ! ! irot(5,1) = 24 ! ! irot(6,1) = 25 ! -! do i = 1,12 -! irot(i,1) = i+6 -! enddo - irot(1,1) = 5 - irot(2,1) = 6 - irot(3,1) = 7 - irot(4,1) = 8 - irot(5,1) = 9 - irot(6,1) = 10 + +!N2 +! irot(1,1) = 5 +! irot(2,1) = 6 +! irot(3,1) = 7 +! irot(4,1) = 8 +! irot(5,1) = 9 +! irot(6,1) = 10 +! +! cmoref(5,1,1) = 1.d0 ! +! cmoref(6,2,1) = 1.d0 ! +! cmoref(7,3,1) = 1.d0 ! +! cmoref(40,4,1) = 1.d0 ! +! cmoref(41,5,1) = 1.d0 ! +! cmoref(42,6,1) = 1.d0 ! +!END N2 + +!HEXATRIENE +! irot(1,1) = 20 +! irot(2,1) = 21 +! irot(3,1) = 22 +! irot(4,1) = 23 +! irot(5,1) = 24 +! irot(6,1) = 25 +! +! cmoref(7,1,1) = 1.d0 ! +! cmoref(26,1,1) = 1.d0 ! +! cmoref(45,2,1) = 1.d0 ! +! cmoref(64,2,1) = 1.d0 ! +! cmoref(83,3,1) = 1.d0 ! +! cmoref(102,3,1) = 1.d0 ! +! cmoref(7,4,1) = 1.d0 ! +! cmoref(26,4,1) = -1.d0 ! +! cmoref(45,5,1) = 1.d0 ! +! cmoref(64,5,1) = -1.d0 ! +! cmoref(83,6,1) = 1.d0 ! +! cmoref(102,6,1) = -1.d0 ! +!END HEXATRIENE + +!!!!H2 H2 CAS +! irot(1,1) = 1 +! irot(2,1) = 2 +! +! cmoref(1,1,1) = 1.d0 +! cmoref(37,2,1) = 1.d0 +!END H2 +!!!! LOCALIZATION ON THE BASIS FUNCTIONS +! do i = 1, nrot(1) +! irot(i,1) = i +! cmoref(i,i,1) = 1.d0 +! enddo + +!END BASISLOC + +! do i = 1, nrot(1) +! irot(i,1) = 4+i +! enddo do i = 1, nrot(1) print*,'irot(i,1) = ',irot(i,1) enddo - pause - cmoref(4,1,1) = 1.d0 ! 2S function - cmoref(5,2,1) = 1.d0 ! 2S function - cmoref(6,3,1) = 1.d0 ! 2S function - cmoref(19,4,1) = 1.d0 ! 2S function - cmoref(20,5,1) = 1.d0 ! 2S function - cmoref(21,6,1) = 1.d0 ! 2S function +! pause ! you define the guess vectors that you want ! the new MO to be close to @@ -138,233 +349,21 @@ ! own guess vectors for the MOs ! The new MOs are provided in output ! in the same order than the guess MOs - - ! C-C bonds - ! 1-2 -! i_atom = 1 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,1,1) = -0.012d0 ! 2S function -! cmoref(2+shift,1,1) = 0.18d0 ! -! cmoref(3+shift,1,1) = 0.1d0 ! - -! cmoref(5+shift,1,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,1,1) = -0.1d0 ! 2pZ function - -! i_atom = 2 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,1,1) = -0.012d0 ! 2S function -! cmoref(2+shift,1,1) = 0.18d0 ! -! cmoref(3+shift,1,1) = 0.1d0 ! - -! cmoref(5+shift,1,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,1,1) = 0.1d0 ! 2pZ function - - -! ! 1-3 -! i_atom = 1 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,2,1) = -0.012d0 ! 2S function -! cmoref(2+shift,2,1) = 0.18d0 ! -! cmoref(3+shift,2,1) = 0.1d0 ! - -! cmoref(5+shift,2,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,2,1) = -0.1d0 ! 2pZ function - -! i_atom = 3 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,2,1) = -0.012d0 ! 2S function -! cmoref(2+shift,2,1) = 0.18d0 ! -! cmoref(3+shift,2,1) = 0.1d0 ! - -! cmoref(5+shift,2,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,2,1) = 0.1d0 ! 2pZ function - -! ! 4-6 -! i_atom = 4 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,3,1) = -0.012d0 ! 2S function -! cmoref(2+shift,3,1) = 0.18d0 ! -! cmoref(3+shift,3,1) = 0.1d0 ! - -! cmoref(5+shift,3,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,3,1) = -0.1d0 ! 2pZ function - -! i_atom = 6 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,3,1) = -0.012d0 ! 2S function -! cmoref(2+shift,3,1) = 0.18d0 ! -! cmoref(3+shift,3,1) = 0.1d0 ! - -! cmoref(5+shift,3,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,3,1) = 0.1d0 ! 2pZ function - - -! ! 6-5 -! i_atom = 6 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,4,1) = -0.012d0 ! 2S function -! cmoref(2+shift,4,1) = 0.18d0 ! -! cmoref(3+shift,4,1) = 0.1d0 ! - -! cmoref(5+shift,4,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,4,1) = 0.1d0 ! 2pZ function - -! i_atom = 5 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,4,1) = -0.012d0 ! 2S function -! cmoref(2+shift,4,1) = 0.18d0 ! -! cmoref(3+shift,4,1) = 0.1d0 ! - -! cmoref(5+shift,4,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,4,1) = -0.1d0 ! 2pZ function - - -! ! 2-4 -! i_atom = 2 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,5,1) = -0.012d0 ! 2S function -! cmoref(2+shift,5,1) = 0.18d0 ! -! cmoref(3+shift,5,1) = 0.1d0 ! - -! cmoref(6+shift,5,1) = 0.1d0 ! 2pZ function - -! i_atom = 4 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,5,1) = -0.012d0 ! 2S function -! cmoref(2+shift,5,1) = 0.18d0 ! -! cmoref(3+shift,5,1) = 0.1d0 ! - -! cmoref(6+shift,5,1) = -0.1d0 ! 2pZ function - - -! ! 3-5 -! i_atom = 3 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,6,1) = -0.012d0 ! 2S function -! cmoref(2+shift,6,1) = 0.18d0 ! -! cmoref(3+shift,6,1) = 0.1d0 ! - -! cmoref(6+shift,6,1) = 0.1d0 ! 2pZ function - -! i_atom = 5 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,6,1) = -0.012d0 ! 2S function -! cmoref(2+shift,6,1) = 0.18d0 ! -! cmoref(3+shift,6,1) = 0.1d0 ! - -! cmoref(6+shift,6,1) = -0.1d0 ! 2pZ function - -! ! C-H bonds -! ! 2-7 -! i_atom = 2 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,7,1) = -0.012d0 ! 2S function -! cmoref(2+shift,7,1) = 0.18d0 ! -! cmoref(3+shift,7,1) = 0.1d0 ! - -! cmoref(5+shift,7,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,7,1) = 0.1d0 ! 2pZ function -! -! i_atom = 7 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,7,1) = 0.12d0 ! 1S function - -! ! 4-10 -! i_atom = 4 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,8,1) = -0.012d0 ! 2S function -! cmoref(2+shift,8,1) = 0.18d0 ! -! cmoref(3+shift,8,1) = 0.1d0 ! - -! cmoref(5+shift,8,1) = -0.1d0 ! 2pX function -! cmoref(6+shift,8,1) = -0.1d0 ! 2pZ function -! -! i_atom = 10 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,8,1) = 0.12d0 ! 1S function - -! ! 5-11 -! i_atom = 5 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,9,1) = -0.012d0 ! 2S function -! cmoref(2+shift,9,1) = 0.18d0 ! -! cmoref(3+shift,9,1) = 0.1d0 ! - -! cmoref(5+shift,9,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,9,1) = -0.1d0 ! 2pZ function -! -! i_atom = 11 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,9,1) = 0.12d0 ! 1S function - -! ! 3-8 -! i_atom = 3 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,10,1) = -0.012d0 ! 2S function -! cmoref(2+shift,10,1) = 0.18d0 ! -! cmoref(3+shift,10,1) = 0.1d0 ! -! -! cmoref(5+shift,10,1) = 0.1d0 ! 2pX function -! cmoref(6+shift,10,1) = 0.1d0 ! 2pZ function -! -! i_atom = 8 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,10,1) = 0.12d0 ! 1S function - -! ! 1-9 -! i_atom = 1 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,11,1) = -0.012d0 ! 2S function -! cmoref(2+shift,11,1) = 0.18d0 ! -! cmoref(3+shift,11,1) = 0.1d0 ! -! -! cmoref(6+shift,11,1) = 0.1d0 ! 2pZ function - -! i_atom = 9 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,11,1) = 0.12d0 ! 1S function - -! -! ! 6-12 -! i_atom = 6 -! shift = (i_atom -1) * 15 -! cmoref(1+shift,12,1) = -0.012d0 ! 2S function -! cmoref(2+shift,12,1) = 0.18d0 ! -! cmoref(3+shift,12,1) = 0.1d0 ! -! -! cmoref(6+shift,12,1) = -0.1d0 ! 2pZ function - -! i_atom = 12 -! shift_h = (6-1) * 15 + (i_atom - 6)*5 -! cmoref(1+shift_h,12,1) = 0.12d0 ! 1S function -! cmoref(12,1,1) = 1.d0 ! - -! cmoref(21,2,1) = 1.d0 ! -! cmoref(30,2,1) = 1.d0 ! - -! cmoref(39,3,1) = 1.d0 ! -! cmoref(48,3,1) = 1.d0 ! - -! cmoref(3,4,1) = 1.d0 ! -! cmoref(12,4,1) =-1.d0 ! - -! cmoref(21,5,1) = 1.d0 ! -! cmoref(30,5,1) =-1.d0 ! - -! cmoref(39,6,1) = 1.d0 ! -! cmoref(48,6,1) =-1.d0 ! +! do i = 1, nrot(1) +! j = 5+(i-1)*15 +! cmoref(j,i,1) = 0.2d0 +! cmoref(j+3,i,1) = 0.12d0 +! print*,'j = ',j +! enddo +! pause print*,'passed the definition of the referent vectors ' - !Building the S (overlap) matrix in the AO basis. - - - do i = 1, ao_num - do j = 1, ao_num - s(i,j,1) = ao_overlap(i,j) + do j =1, ao_num + s(i,j,1) = ao_overlap(i,j) enddo enddo !Now big loop over symmetry @@ -398,20 +397,13 @@ ! do i=1,nmo(isym) - do i=1,ao_num - do j=1,nrot(isym) - - ddum(i,j)=0.d0 - - do k=1,ao_num - - ddum(i,j)=ddum(i,j)+s(i,k,isym)*cmo(k,irot(j,isym),isym) - - enddo - - enddo - + do i=1,ao_num + ddum(i,j)=0.d0 + do k=1,ao_num + ddum(i,j)=ddum(i,j)+s(i,k,isym)*cmo(k,irot(j,isym),isym) + enddo + enddo enddo @@ -441,7 +433,7 @@ do i=1,nrot(isym) do j=1,ao_num - write (6,*) 'isym,',isym,nrot(isym),nmo(isym) +! write (6,*) 'isym,',isym,nrot(isym),nmo(isym) newcmo(j,irot(i,isym),isym)=0.d0 do k=1,nrot(isym) newcmo(j,irot(i,isym),isym)=newcmo(j,irot(i,isym),isym) + cmo(j,irot(k,isym),isym)*t(k,i) @@ -459,7 +451,7 @@ enddo !big loop over symmetry - 10 format (4E20.12) + 10 format (4E18.12) ! Now we copyt the newcmo into the mo_coef @@ -472,9 +464,7 @@ enddo enddo enddo -! if(dabs(newcmo(3,19,1) - mo_coef(3,19)) .gt.1.d-10 )then - print*,'mo_coef(3,19)',mo_coef(3,19) - pause +! pause ! we say that it hase been touched, and valid and that everything that diff --git a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py index a140ec1e..e911af28 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py @@ -264,7 +264,7 @@ def print_mo_coef(mo_coef_block, l_l_sym): i_a = int(l[1]) - 1 sym = l[2] - print l_label[i_a], sym, " ".join('{: 3.8f}'.format(i) + print l_label[i_a], sym, " ".join('{0: 3.8f}'.format(i) for i in a[i]) if i_block != nb_block - 1: diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index e5006b12..ae0064cf 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -8,11 +8,22 @@ copy_buffer declarations decls_main deinit_thread -do_double_excitations +skip +init_main +filter_integrals +filter2p +filter2h2p_double +filter2h2p_single filter1h filter1p -filter2h2p -filter2p +only_2p_single +only_2p_double +filter_only_1h1p_single +filter_only_1h1p_double +filter_only_1h2p_single +filter_only_1h2p_double +filter_only_2h2p_single +filter_only_2h2p_double filterhole filter_integrals filter_only_1h1p_double @@ -182,7 +193,7 @@ class H_apply(object): if (is_a_2p(hole)) cycle """ def filter_1p(self): - self["filter0p"] = """ + self["filter1p"] = """ ! ! DIR$ FORCEINLINE if (is_a_1p(hole)) cycle """ @@ -208,6 +219,27 @@ class H_apply(object): if (is_a_1h1p(key).eqv..False.) cycle """ + def filter_only_2h2p(self): + self["filter_only_2h2p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_two_holes_two_particles(hole).eqv..False.) cycle + """ + self["filter_only_1h1p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_two_holes_two_particles(key).eqv..False.) cycle + """ + + + def filter_only_1h2p(self): + self["filter_only_1h2p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h2p(hole).eqv..False.) cycle + """ + self["filter_only_1h2p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h2p(key).eqv..False.) cycle + """ + def unset_skip(self): self["skip"] = """ @@ -215,9 +247,12 @@ class H_apply(object): def set_filter_2h_2p(self): - self["filter2h2p"] = """ + self["filter2h2p_double"] = """ if (is_a_two_holes_two_particles(key)) cycle """ + self["filter2h2p_single"] = """ + if (is_a_two_holes_two_particles(hole)) cycle + """ def set_perturbation(self,pert): diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 67ccf03c..8d420b15 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -206,3 +206,176 @@ BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)] l_to_charater(4)='G' END_PROVIDER + + BEGIN_PROVIDER [ integer, Nucl_N_Aos, (nucl_num)] +&BEGIN_PROVIDER [ integer, N_AOs_max ] + implicit none + integer :: i + BEGIN_DOC + ! Number of AOs per atom + END_DOC + Nucl_N_Aos = 0 + do i = 1, ao_num + Nucl_N_Aos(ao_nucl(i)) +=1 + enddo + N_AOs_max = maxval(Nucl_N_Aos) +END_PROVIDER + + BEGIN_PROVIDER [ integer, Nucl_Aos, (nucl_num,N_AOs_max)] + implicit none + BEGIN_DOC + ! List of AOs attached on each atom + END_DOC + integer :: i + integer, allocatable :: nucl_tmp(:) + allocate(nucl_tmp(nucl_num)) + nucl_tmp = 0 + Nucl_Aos = 0 + do i = 1, ao_num + nucl_tmp(ao_nucl(i))+=1 + Nucl_Aos(ao_nucl(i),nucl_tmp(ao_nucl(i))) = i + enddo + deallocate(nucl_tmp) +END_PROVIDER + + + BEGIN_PROVIDER [ integer, Nucl_list_shell_Aos, (nucl_num,N_AOs_max)] +&BEGIN_PROVIDER [ integer, Nucl_num_shell_Aos, (nucl_num)] + implicit none + integer :: i,j,k + BEGIN_DOC + ! Index of the shell type Aos and of the corresponding Aos + ! Per convention, for P,D,F and G AOs, we take the index + ! of the AO with the the corresponding power in the "X" axis + END_DOC + do i = 1, nucl_num + Nucl_num_shell_Aos(i) = 0 + + do j = 1, Nucl_N_Aos(i) + if(ao_l(Nucl_Aos(i,j))==0)then + ! S type function + Nucl_num_shell_Aos(i)+=1 + Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) + elseif(ao_l(Nucl_Aos(i,j))==1)then + ! P type function + if(ao_power(Nucl_Aos(i,j),1)==1)then + Nucl_num_shell_Aos(i)+=1 + Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) + endif + elseif(ao_l(Nucl_Aos(i,j))==2)then + ! D type function + if(ao_power(Nucl_Aos(i,j),1)==2)then + Nucl_num_shell_Aos(i)+=1 + Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) + endif + elseif(ao_l(Nucl_Aos(i,j))==3)then + ! F type function + if(ao_power(Nucl_Aos(i,j),1)==3)then + Nucl_num_shell_Aos(i)+=1 + Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) + endif + elseif(ao_l(Nucl_Aos(i,j))==4)then + ! G type function + if(ao_power(Nucl_Aos(i,j),1)==4)then + Nucl_num_shell_Aos(i)+=1 + Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) + endif + endif + + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] + implicit none + integer :: i + character*(4) :: give_ao_character_space + do i=1,ao_num + + if(ao_l(i)==0)then + ! S type AO + give_ao_character_space = 'S ' + elseif(ao_l(i) == 1)then + ! P type AO + if(ao_power(i,1)==1)then + give_ao_character_space = 'X ' + elseif(ao_power(i,2) == 1)then + give_ao_character_space = 'Y ' + else + give_ao_character_space = 'Z ' + endif + elseif(ao_l(i) == 2)then + ! D type AO + if(ao_power(i,1)==2)then + give_ao_character_space = 'XX ' + elseif(ao_power(i,2) == 2)then + give_ao_character_space = 'YY ' + elseif(ao_power(i,3) == 2)then + give_ao_character_space = 'ZZ ' + elseif(ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'XY ' + elseif(ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'XZ ' + else + give_ao_character_space = 'YZ ' + endif + elseif(ao_l(i) == 3)then + ! F type AO + if(ao_power(i,1)==3)then + give_ao_character_space = 'XXX ' + elseif(ao_power(i,2) == 3)then + give_ao_character_space = 'YYY ' + elseif(ao_power(i,3) == 3)then + give_ao_character_space = 'ZZZ ' + elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'XXY ' + elseif(ao_power(i,1) == 2 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'XXZ ' + elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1)then + give_ao_character_space = 'YYX ' + elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'YYZ ' + elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1)then + give_ao_character_space = 'ZZX ' + elseif(ao_power(i,3) == 2 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'ZZY ' + elseif(ao_power(i,3) == 1 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'XYZ ' + endif + elseif(ao_l(i) == 4)then + ! G type AO + if(ao_power(i,1)==4)then + give_ao_character_space = 'XXXX' + elseif(ao_power(i,2) == 4)then + give_ao_character_space = 'YYYY' + elseif(ao_power(i,3) == 4)then + give_ao_character_space = 'ZZZZ' + elseif(ao_power(i,1) == 3 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'XXXY' + elseif(ao_power(i,1) == 3 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'XXXZ' + elseif(ao_power(i,2) == 3 .and. ao_power(i,1) == 1)then + give_ao_character_space = 'YYYX' + elseif(ao_power(i,2) == 3 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'YYYZ' + elseif(ao_power(i,3) == 3 .and. ao_power(i,1) == 1)then + give_ao_character_space = 'ZZZX' + elseif(ao_power(i,3) == 3 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'ZZZY' + elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 2)then + give_ao_character_space = 'XXYY' + elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 2)then + give_ao_character_space = 'YYZZ' + elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'XXYZ' + elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then + give_ao_character_space = 'YYXZ' + elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then + give_ao_character_space = 'ZZXY' + endif + endif + ao_l_char_space(i) = give_ao_character_space + enddo +END_PROVIDER diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 4441fb22..4984d9a8 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -212,6 +212,12 @@ logical function is_a_two_holes_two_particles(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: i,i_diff + integer :: number_of_holes, number_of_particles + is_a_two_holes_two_particles = .False. + if(number_of_holes(key_in) == 2 .and. number_of_particles(key_in) == 2)then + is_a_two_holes_two_particles = .True. + return + endif i_diff = 0 if(N_int == 1)then i_diff = i_diff & @@ -456,6 +462,17 @@ logical function is_a_1h1p(key_in) end +logical function is_a_1h2p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h2p = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.2)then + is_a_1h2p = .True. + endif + +end + logical function is_a_1h(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 6fe36c57..7bb6e16e 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -95,9 +95,40 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ] END_PROVIDER +BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ] + implicit none + BEGIN_DOC + ! Number of bitmasks for generators + END_DOC + logical :: exists + PROVIDE ezfio_filename + + call ezfio_has_bitmasks_N_mask_gen(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_generators_bitmask_restart = 1 + endif + ASSERT (N_generators_bitmask_restart > 0) + +END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask) ] + + +BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ] implicit none BEGIN_DOC ! Bitmasks for generator determinants. @@ -306,7 +337,7 @@ END_PROVIDER n_inact_orb = 0 n_virt_orb = 0 - if(N_generators_bitmask == 1)then + if(N_generators_bitmask_restart == 1)then do j = 1, N_int inact_bitmask(j,1) = xor(generators_bitmask_restart(j,1,1,1),cas_bitmask(j,1,1)) inact_bitmask(j,2) = xor(generators_bitmask_restart(j,2,1,1),cas_bitmask(j,2,1)) @@ -319,15 +350,15 @@ END_PROVIDER i_hole = 1 i_gen = 1 do i = 1, N_int - inact_bitmask(i,1) = generators_bitmask(i,1,i_hole,i_gen) - inact_bitmask(i,2) = generators_bitmask(i,2,i_hole,i_gen) + inact_bitmask(i,1) = generators_bitmask_restart(i,1,i_hole,i_gen) + inact_bitmask(i,2) = generators_bitmask_restart(i,2,i_hole,i_gen) n_inact_orb += popcnt(inact_bitmask(i,1)) enddo i_part = 2 i_gen = 3 do i = 1, N_int - virt_bitmask(i,1) = generators_bitmask(i,1,i_part,i_gen) - virt_bitmask(i,2) = generators_bitmask(i,2,i_part,i_gen) + virt_bitmask(i,1) = generators_bitmask_restart(i,1,i_part,i_gen) + virt_bitmask(i,2) = generators_bitmask_restart(i,2,i_part,i_gen) n_virt_orb += popcnt(virt_bitmask(i,1)) enddo endif diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 45f08f75..5d9198c4 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -170,6 +170,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl logical :: check_double_excitation logical :: is_a_1h1p + logical :: is_a_1h2p logical :: is_a_1h logical :: is_a_1p logical :: is_a_2p @@ -299,8 +300,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) - $filter2h2p + $filter2h2p_double $filter_only_1h1p_double + $filter_only_1h2p_double + $filter_only_2h2p_double $only_2p_double key_idx += 1 do k=1,N_int @@ -348,8 +351,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,ispin) = ibset(key(k,ispin),l) - $filter2h2p + $filter2h2p_double $filter_only_1h1p_double + $filter_only_1h2p_double + $filter_only_2h2p_double $only_2p_double key_idx += 1 do k=1,N_int @@ -418,6 +423,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato logical :: check_double_excitation logical :: is_a_1h1p + logical :: is_a_1h2p logical :: is_a_1h logical :: is_a_1p logical :: is_a_2p @@ -494,8 +500,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato $filter1h $filter1p $filter2p - $filter2h2p + $filter2h2p_single $filter_only_1h1p_single + $filter_only_1h2p_single + $filter_only_2h2p_single key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) @@ -521,4 +529,3 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato end - diff --git a/src/Determinants/SC2.irp.f b/src/Determinants/SC2.irp.f index ea942307..4f321b87 100644 --- a/src/Determinants/SC2.irp.f +++ b/src/Determinants/SC2.irp.f @@ -1,4 +1,4 @@ -subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) +subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence) use bitmasks implicit none BEGIN_DOC @@ -21,6 +21,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(inout) :: u_in(dim_in,N_st) double precision, intent(out) :: energies(N_st) + double precision, intent(out) :: diag_H_elements(dim_in) double precision, intent(in) :: convergence ASSERT (N_st > 0) ASSERT (sze > 0) @@ -197,6 +198,9 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) converged = dabs(e_corr_double - e_corr_double_before) < convergence converged = converged if (converged) then + do i = 1, dim_in + diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i) + enddo exit endif e_corr_double_before = e_corr_double diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 4476ed45..a7727cda 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -58,7 +58,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ] else psi_det_size = 1 endif - psi_det_size = max(psi_det_size,10000) + psi_det_size = max(psi_det_size,100000) call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays') END_PROVIDER diff --git a/src/Determinants/diagonalize_CI_SC2.irp.f b/src/Determinants/diagonalize_CI_SC2.irp.f index 97161ad3..498792d9 100644 --- a/src/Determinants/diagonalize_CI_SC2.irp.f +++ b/src/Determinants/diagonalize_CI_SC2.irp.f @@ -23,8 +23,10 @@ END_PROVIDER threshold_convergence_SC2 = 1.d-10 END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ] implicit none BEGIN_DOC ! Eigenvectors/values of the CI matrix @@ -39,7 +41,8 @@ END_PROVIDER enddo call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & - size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) +! size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) END_PROVIDER subroutine diagonalize_CI_SC2 @@ -54,5 +57,6 @@ subroutine diagonalize_CI_SC2 psi_coef(i,j) = CI_SC2_eigenvectors(i,j) enddo enddo - SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors + SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2 +! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors end diff --git a/src/Determinants/save_natorb.irp.f b/src/Determinants/save_natorb.irp.f index e56f9821..674ba32e 100644 --- a/src/Determinants/save_natorb.irp.f +++ b/src/Determinants/save_natorb.irp.f @@ -2,5 +2,6 @@ program save_natorb read_wf = .True. touch read_wf call save_natural_mos + call save_ref_determinant end diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 020b5dfa..fdcf4038 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -230,7 +230,6 @@ subroutine clear_ao_map end - !! MO Map !! ====== diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 630ef699..0a468c24 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -70,7 +70,7 @@ subroutine add_integrals_to_map(mask_ijkl) integer :: i2,i3,i4 double precision,parameter :: thr_coef = 1.d-10 - PROVIDE ao_bielec_integrals_in_map + PROVIDE ao_bielec_integrals_in_map mo_coef !Get list of MOs for i,j,k and l !------------------------------- @@ -328,7 +328,7 @@ end double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) if (.not.do_direct_integrals) then - PROVIDE ao_bielec_integrals_in_map + PROVIDE ao_bielec_integrals_in_map mo_coef endif mo_bielec_integral_jj_from_ao = 0.d0 @@ -494,4 +494,13 @@ subroutine clear_mo_map call map_deinit(mo_integrals_map) FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + + +end + +subroutine provide_all_mo_integrals + implicit none + provide mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti + provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + end diff --git a/src/Integrals_Monoelec/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f index 714222ec..5bae9868 100644 --- a/src/Integrals_Monoelec/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -5,6 +5,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to ! array of the mono electronic hamiltonian on the MOs basis ! : sum of the kinetic and nuclear electronic potential END_DOC + print*,'Providing the mono electronic integrals' do j = 1, mo_tot_num do i = 1, mo_tot_num mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)