From 3b8239976ab96929f652dfc363dcfc5d62b44a64 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 20 Mar 2017 16:21:00 +0100 Subject: [PATCH] minor modifs in FOBOCI --- plugins/DDCI_selected/NEEDED_CHILDREN_MODULES | 2 +- plugins/DDCI_selected/ddci.irp.f | 2 +- plugins/FOBOCI/density_matrix.irp.f | 46 ++++++++------ plugins/FOBOCI/dress_simple.irp.f | 62 ++++++++++--------- .../foboci_lmct_mlct_threshold_old.irp.f | 23 +++---- plugins/FOBOCI/generators_restart_save.irp.f | 20 ++++-- plugins/FOBOCI/routines_foboci.irp.f | 48 +++++++------- plugins/MRPT/MRPT_Utils.main.irp.f | 8 +-- plugins/MRPT_Utils/energies_cas.irp.f | 32 ---------- plugins/MRPT_Utils/mrpt_utils.irp.f | 9 +-- ...gonalize_restart_and_save_all_states.irp.f | 2 +- ...gonalize_restart_and_save_two_states.irp.f | 27 -------- 12 files changed, 121 insertions(+), 160 deletions(-) delete mode 100644 src/Determinants/diagonalize_restart_and_save_two_states.irp.f diff --git a/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES b/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES index 0b7ce8a9..d212e150 100644 --- a/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES +++ b/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_CAS Davidson +Perturbation Selectors_full Generators_CAS Davidson Psiref_CAS diff --git a/plugins/DDCI_selected/ddci.irp.f b/plugins/DDCI_selected/ddci.irp.f index 0bfb324f..a1824857 100644 --- a/plugins/DDCI_selected/ddci.irp.f +++ b/plugins/DDCI_selected/ddci.irp.f @@ -5,7 +5,7 @@ program ddci double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) integer :: N_st, degree - N_st = N_states_diag + N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) character*(64) :: perturbation diff --git a/plugins/FOBOCI/density_matrix.irp.f b/plugins/FOBOCI/density_matrix.irp.f index 5a06d5d7..14a2fefa 100644 --- a/plugins/FOBOCI/density_matrix.irp.f +++ b/plugins/FOBOCI/density_matrix.irp.f @@ -1,12 +1,12 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_generators_restart, (mo_tot_num_align,mo_tot_num) ] &BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_generators_restart, (mo_tot_num_align,mo_tot_num) ] -&BEGIN_PROVIDER [ double precision, norm_generators_restart] +&BEGIN_PROVIDER [ double precision, norm_generators_restart, (N_states)] implicit none BEGIN_DOC ! Alpha and beta one-body density matrix for the generators restart END_DOC - integer :: j,k,l,m + integer :: j,k,l,m,istate integer :: occ(N_int*bit_kind_size,2) double precision :: ck, cl, ckl double precision :: phase @@ -14,29 +14,37 @@ integer :: exc(0:2,2,2),n_occ_alpha double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) integer :: degree_respect_to_HF_k - integer :: degree_respect_to_HF_l,index_ref_generators_restart - double precision :: inv_coef_ref_generators_restart + integer :: degree_respect_to_HF_l,index_ref_generators_restart(N_states) + double precision :: inv_coef_ref_generators_restart(N_states) integer :: i print*, 'providing the one_body_dm_mo_alpha_generators_restart' - do i = 1, N_det_generators_restart - ! Find the reference determinant for intermediate normalization - call get_excitation_degree(ref_generators_restart,psi_det_generators_restart(1,1,i),degree,N_int) - if(degree == 0)then - index_ref_generators_restart = i - inv_coef_ref_generators_restart = 1.d0/psi_coef_generators_restart(i,1) - exit - endif + do istate = 1, N_states + do i = 1, N_det_generators_restart + ! Find the reference determinant for intermediate normalization + call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det_generators_restart(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart(istate) = i + inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef_generators_restart(i,istate) + exit + endif + enddo enddo norm_generators_restart = 0.d0 - do i = 1, N_det_generators_restart - psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart - norm_generators_restart += psi_coef_generators_restart(i,1)**2 + do istate = 1, N_states + do i = 1, N_det_generators_restart + psi_coef_generators_restart(i,istate) = psi_coef_generators_restart(i,istate) * inv_coef_ref_generators_restart(istate) + norm_generators_restart(istate) += psi_coef_generators_restart(i,istate)**2 + enddo enddo - double precision :: inv_norm - inv_norm = 1.d0/dsqrt(norm_generators_restart) - do i = 1, N_det_generators_restart - psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_norm + double precision :: inv_norm(N_States) + do istate = 1, N_states + inv_norm(istate) = 1.d0/dsqrt(norm_generators_restart(istate)) + enddo + do istate = 1, N_states + do i = 1, N_det_generators_restart + psi_coef_generators_restart(i,istate) = psi_coef_generators_restart(i,istate) * inv_norm(istate) + enddo enddo diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index fabbd834..c74d08e7 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -107,7 +107,6 @@ subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_per !enddo !soft_touch psi_selectors psi_selectors_coef !if(do_it_perturbative)then - print*, 'is_ok_perturbative',is_ok_perturbative if(is_ok.or.is_ok_perturbative)then N_det = N_det_generators do m = 1, N_states @@ -117,7 +116,6 @@ subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_per psi_det(l,2,k) = psi_det_generators_input(l,2,k) enddo psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m) - print*, 'psi_coef(k,m)',psi_coef(k,m) enddo enddo soft_touch psi_det psi_coef N_det @@ -150,7 +148,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators) - integer :: i,j,degree,index_ref_generators_restart,i_count,k,i_det_no_ref + integer :: i,j,degree,index_ref_generators_restart(N_states),i_count,k,i_det_no_ref double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average logical :: is_a_ref_det(Ndet_generators) @@ -168,17 +166,19 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener enddo + integer :: istate + do istate = 1, N_states + do i = 1, Ndet_generators + call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det_generators_input(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart(istate) = i + exit + endif + enddo + enddo do i = 1, Ndet_generators - call get_excitation_degree(ref_generators_restart,psi_det_generators_input(1,1,i),degree,N_int) - if(degree == 0)then - index_ref_generators_restart = i - endif do j = 1, Ndet_generators call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij) ! Fill the zeroth order H matrix - if(i==j)then - call debug_det(psi_det_generators_input(1,1,i),N_int) - print*, hij - endif dressed_H_matrix(i,j) = hij enddo enddo @@ -189,15 +189,21 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener i_det_no_ref +=1 diag_h_mat_average+=dressed_H_matrix(i,i) enddo + double precision :: average_ref_h_mat + average_ref_h_mat = 0.d0 + do istate = 1, N_states + average_ref_h_mat += dressed_H_matrix(index_ref_generators_restart(istate),index_ref_generators_restart(istate)) + enddo + average_ref_h_mat = 1.d0/dble(N_states) diag_h_mat_average = diag_h_mat_average/dble(i_det_no_ref) print*,'diag_h_mat_average = ',diag_h_mat_average - print*,'ref h_mat = ',dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) + print*,'ref h_mat average = ',average_ref_h_mat integer :: number_of_particles, number_of_holes ! Filter the the MLCT that are higher than 27.2 eV in energy with respect to the reference determinant do i = 1, Ndet_generators if(is_a_ref_det(i))cycle if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then - if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then + if(diag_h_mat_average - average_ref_h_mat .gt.2.d0)then is_ok = .False. exit_loop = .True. return @@ -206,7 +212,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener ! Filter the the LMCT that are higher than 54.4 eV in energy with respect to the reference determinant if(number_of_holes(psi_det_generators_input(1,1,i)).eq.1 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.0)then - if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then + if(diag_h_mat_average - average_ref_h_mat .gt.1.d0)then is_ok = .False. return endif @@ -214,7 +220,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener exit enddo - call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix + call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the naked matrix double precision :: s2(N_det_generators),E_ref(N_states) integer :: i_state(N_states) @@ -238,18 +244,12 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener do i = 1, N_states i_state(i) = i E_ref(i) = eigvalues(i) - print*, 'E_ref(i)',E_ref(i) enddo endif - do i = 1,N_states - print*,'i_state = ',i_state(i) - enddo do k = 1, N_states - print*,'state ',k do i = 1, Ndet_generators - psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k)) + psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart(k),i_state(k)) psi_coef_ref(i,k) = eigvectors(i,i_state(k)) - print*,'psi_coef_ref(i) = ',psi_coef_ref(i,k) enddo enddo if(verbose)then @@ -262,7 +262,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener do k = 1, N_states print*,'state ',k do i = 1, Ndet_generators - print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i) + print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart(k),index_ref_generators_restart(k)),is_a_ref_det(i) enddo enddo endif @@ -283,18 +283,20 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix integer :: i_good_state(0:N_states) i_good_state(0) = 0 - do i = 1, Ndet_generators + do k = 1, N_states +! print*,'state',k + do i = 1, Ndet_generators ! State following - do k = 1, N_states accu = 0.d0 do j =1, Ndet_generators - print*,'',eigvectors(j,i) , psi_coef_ref(j,k) accu += eigvectors(j,i) * psi_coef_ref(j,k) enddo - print*,'accu = ',accu +! print*,i,accu if(dabs(accu).ge.0.60d0)then i_good_state(0) +=1 i_good_state(i_good_state(0)) = i + print*, 'state, ovrlap',k,i,accu + exit endif enddo if(i_good_state(0)==N_states)then @@ -309,14 +311,14 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener accu = 0.d0 do k = 1, N_states do i = 1, Ndet_generators - psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k)) + psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart(k),i_state(k)) enddo enddo if(verbose)then do k = 1, N_states print*,'state ',k do i = 1, Ndet_generators - print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i) + print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart(k),index_ref_generators_restart(k)),is_a_ref_det(i) enddo enddo endif @@ -338,7 +340,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener do i = 1, Ndet_generators if(is_a_ref_det(i))cycle do k = 1, N_states - print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative +! print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold_perturbative)then is_ok_perturbative = .False. exit diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index 5b549fc3..746704c2 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -46,6 +46,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter) lmct = .True. integer :: i_hole_osoci i_hole_osoci = list_inact(i) +! if(i_hole_osoci.ne.26)cycle print*,'--------------------------' ! First set the current generators to the one of restart call check_symetry(i_hole_osoci,thr,test_sym) @@ -55,11 +56,6 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter) print*,'i_hole_osoci = ',i_hole_osoci call create_restart_and_1h(i_hole_osoci) call set_generators_to_psi_det - print*,'Passed set generators' - integer :: m - do m = 1, N_det_generators - call debug_det(psi_det_generators(1,1,m),N_int) - enddo call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) double precision :: e_pt2 @@ -88,9 +84,9 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter) call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single(e_pt2) ! call make_s2_eigenfunction_first_order - threshold_davidson = 1.d-6 - soft_touch threshold_davidson davidson_criterion - call diagonalize_ci +! threshold_davidson = 1.d-6 +! soft_touch threshold_davidson davidson_criterion +! call diagonalize_ci double precision :: hkl call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) hkl = dressing_matrix(1,1) @@ -123,6 +119,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter) do i = 1, n_virt_orb integer :: i_particl_osoci i_particl_osoci = list_virt(i) +! cycle print*,'--------------------------' ! First set the current generators to the one of restart @@ -157,11 +154,11 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter) enddo enddo call all_single(e_pt2) - call make_s2_eigenfunction_first_order - threshold_davidson = 1.d-6 - soft_touch threshold_davidson davidson_criterion - - call diagonalize_ci +! call make_s2_eigenfunction_first_order +! threshold_davidson = 1.d-6 +! soft_touch threshold_davidson davidson_criterion +! +! call diagonalize_ci deallocate(dressing_matrix) else if(exit_loop)then diff --git a/plugins/FOBOCI/generators_restart_save.irp.f b/plugins/FOBOCI/generators_restart_save.irp.f index 669c899d..6ec528cf 100644 --- a/plugins/FOBOCI/generators_restart_save.irp.f +++ b/plugins/FOBOCI/generators_restart_save.irp.f @@ -21,23 +21,19 @@ END_PROVIDER 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 [ integer(bit_kind), ref_generators_restart, (N_int,2,N_states) ] &BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ] implicit none BEGIN_DOC ! read wf ! END_DOC - integer :: i, k + integer :: i, k,j integer, save :: ifirst = 0 double precision, allocatable :: psi_coef_read(:,:) print*, ' Providing psi_det_generators_restart' 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 @@ -45,6 +41,18 @@ END_PROVIDER psi_coef_generators_restart(i,k) = psi_coef_read(i,k) enddo enddo + do k = 1, N_states + do i = 1, N_det_generators_restart + if(dabs(psi_coef_generators_restart(i,k)).gt.0.5d0)then + do j = 1, N_int + ref_generators_restart(j,1,k) = psi_det_generators_restart(j,1,i) + ref_generators_restart(j,2,k) = psi_det_generators_restart(j,2,i) + enddo + exit + endif + enddo + call debug_det(ref_generators_restart(1,1,k),N_int) + enddo ifirst = 1 deallocate(psi_coef_read) else diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index cda7dd75..db683c96 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -2,7 +2,7 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole) implicit none integer, intent(in) :: i_hole double precision, intent(out) :: norm(N_states) - integer :: i,j,degree,index_ref_generators_restart,k + integer :: i,j,degree,index_ref_generators_restart(N_states),k integer:: number_of_holes,n_h, number_of_particles,n_p integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:) integer, allocatable :: index_one_p(:) @@ -24,16 +24,18 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole) n_good_hole = 0 ! Find the one holes and one hole one particle is_a_ref_det = .False. + integer :: istate + do istate = 1, N_States + do i = 1, N_det + ! Find the reference determinant for intermediate normalization + call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart(istate) = i + inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef(i,istate) + endif + enddo + enddo do i = 1, N_det - ! Find the reference determinant for intermediate normalization - call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int) - if(degree == 0)then - index_ref_generators_restart = i - do k = 1, N_states - inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k) - enddo - endif - ! Find all the determinants present in the reference wave function do j = 1, N_det_generators_restart call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int) @@ -67,7 +69,7 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole) do k = 1,N_states print*,'state ',k do i = 1, n_good_hole - print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart,k) + print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart(k),k) enddo print*,'' enddo @@ -110,7 +112,7 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl) implicit none integer, intent(in) :: i_particl double precision, intent(out) :: norm(N_states) - integer :: i,j,degree,index_ref_generators_restart,k + integer :: i,j,degree,index_ref_generators_restart(N_states),k integer:: number_of_holes,n_h, number_of_particles,n_p integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:) integer, allocatable :: index_one_p(:),index_one_hole_two_p(:) @@ -139,16 +141,18 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl) ! Find the one holes and one hole one particle i_count = 0 is_a_ref_det = .False. - do i = 1, N_det - call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int) - if(degree == 0)then - index_ref_generators_restart = i - do k = 1, N_states - inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k) - enddo -! cycle - endif + integer :: istate + do istate = 1, N_states + do i = 1, N_det + call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart(istate) = i + inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef(i,istate) + endif + enddo + enddo + do i = 1, N_det ! Find all the determinants present in the reference wave function do j = 1, N_det_generators_restart call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int) @@ -184,7 +188,7 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl) do k = 1, N_states print*,'state ',k do i = 1, n_good_particl - print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart,k) + print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart(k),k) enddo print*,'' enddo diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f index d78fd724..3b559232 100644 --- a/plugins/MRPT/MRPT_Utils.main.irp.f +++ b/plugins/MRPT/MRPT_Utils.main.irp.f @@ -32,13 +32,13 @@ subroutine routine_3 call save_wavefunction_general(N_det_ref,N_states_diag_heff,psi_ref,N_det_ref,CI_dressed_pt2_new_eigenvectors) endif if(N_states.gt.1)then - print*, 'Energy differences : E(0) - E(i)' + print*, 'Energy differences : E(i) - E(0)' do i = 2, N_States print*,'State',i write(*,'(A12,X,I3,A3,XX,F20.16)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i) - write(*,'(A12,X,I3,A3,XX,F20.16)') 'Variational ', i,' = ', psi_ref_average_value(1) - psi_ref_average_value(i) - write(*,'(A12,X,I3,A3,XX,F20.16)') 'Perturbative', i,' = ', psi_ref_average_value(1)+second_order_pt_new(1) - (psi_ref_average_value(i)+second_order_pt_new(i)) - write(*,'(A12,X,I3,A3,XX,F20.16)') 'Dressed ', i,' = ', CI_dressed_pt2_new_energy(1) - CI_dressed_pt2_new_energy(i) + write(*,'(A12,X,I3,A3,XX,F20.16)') 'Variational ', i,' = ', -(psi_ref_average_value(1) - psi_ref_average_value(i)) + write(*,'(A12,X,I3,A3,XX,F20.16)') 'Perturbative', i,' = ', -(psi_ref_average_value(1)+second_order_pt_new(1) - (psi_ref_average_value(i)+second_order_pt_new(i))) + write(*,'(A12,X,I3,A3,XX,F20.16)') 'Dressed ', i,' = ', -( CI_dressed_pt2_new_energy(1) - CI_dressed_pt2_new_energy(i) ) enddo endif diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index f7e48e4f..02077ebc 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -688,10 +688,6 @@ END_PROVIDER call get_mono_excitation(psi_in_out(1,1,i),psi_ref(1,1,i),exc,phase,N_int) do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j)* hij * phase -! if(orb_i == 5 .and. orb_v == 20)then - if(orb_i == 2 .and. orb_v == 6 )then - print*, i, ispin - endif norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo enddo @@ -702,12 +698,7 @@ END_PROVIDER one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0 else norm_no_inv(j,ispin) = norm(j,ispin) -! one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin) norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) -! if(orb_i == 5 .and. orb_v == 20)then - if(orb_i == 2 .and. orb_v == 6 )then - print*,ispin ,norm(j,ispin) - endif endif enddo integer :: iorb_annil,hole_particle,spin_exc,orb @@ -721,12 +712,8 @@ END_PROVIDER do i = 1, N_det_ref do j = 1, N_int - ! psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1)) - ! psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1)) psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,2,i) = psi_active(j,2,i) - ! psi_in_out(j,1,i) = psi_ref(j,1,i) - ! psi_in_out(j,2,i) = psi_ref(j,2,i) enddo enddo do state_target = 1, N_states @@ -740,30 +727,11 @@ END_PROVIDER do state_target = 1, N_states if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - & -! 0.5d0 * (energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2)) ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) else one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0 endif - if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then -! if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.29269686324d0) < 1.d-11)then - print*, '' - print*, orb_i,orb_v - print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) !/ norm_bis(state_target,1) - print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) !/ norm_bis(state_target,2) - print*, fock_core_inactive_total_spin_trace(orb_i,1) - print*, fock_virt_total_spin_trace(orb_v,1) - print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target) - print*, '' - endif - if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .gt. 1.d-10)then - write(*,'(F11.8)'), one_anhil_one_creat_inact_virt(iorb,vorb,state_target) -! if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .lt. 1.d-2)then -! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0 -! print*, orb_i,orb_v -! endif - endif enddo enddo enddo diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 5cf67eb6..52bd258b 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -305,8 +305,8 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag_heff) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag_heff) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag_heff) ] + &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states) ] + &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states) ] BEGIN_DOC ! Eigenvectors/values of the CI matrix END_DOC @@ -394,6 +394,7 @@ END_PROVIDER psi_tmp(i) = eigenvectors(i,iorder(1)) enddo CI_electronic_dressed_pt2_new_energy(i_state) = eigenvalues(iorder(1)) + print*, 'CI_electronic_dressed_pt2_new_energy',CI_electronic_dressed_pt2_new_energy(i_state) call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det_ref,psi_det,N_int,1,N_det_ref) print*,'S^2 = ', CI_dressed_pt2_new_eigenvectors_s2(i_state) enddo @@ -502,7 +503,7 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag_heff) ] +BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states) ] implicit none BEGIN_DOC ! N_states lowest eigenvalues of the CI matrix @@ -511,7 +512,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag_hef integer :: j character*(8) :: st call write_time(output_determinants) - do j=1,N_states_diag_heff + do j=1,N_states CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st)) diff --git a/src/Davidson/diagonalize_restart_and_save_all_states.irp.f b/src/Davidson/diagonalize_restart_and_save_all_states.irp.f index 3bdc37c5..393ff63a 100644 --- a/src/Davidson/diagonalize_restart_and_save_all_states.irp.f +++ b/src/Davidson/diagonalize_restart_and_save_all_states.irp.f @@ -9,7 +9,7 @@ subroutine routine implicit none call diagonalize_CI print*,'N_det = ',N_det - call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) + call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) diff --git a/src/Determinants/diagonalize_restart_and_save_two_states.irp.f b/src/Determinants/diagonalize_restart_and_save_two_states.irp.f deleted file mode 100644 index 97fed531..00000000 --- a/src/Determinants/diagonalize_restart_and_save_two_states.irp.f +++ /dev/null @@ -1,27 +0,0 @@ -program diag_and_save - implicit none - read_wf = .True. - touch read_wf - call routine -end - -subroutine routine - implicit none - integer :: igood_state_1,igood_state_2 - double precision, allocatable :: psi_coef_tmp(:,:) - integer :: i - print*,'N_det = ',N_det -!call diagonalize_CI - write(*,*)'Which couple of states would you like to save ?' - read(5,*)igood_state_1,igood_state_2 - allocate(psi_coef_tmp(n_det,2)) - do i = 1, N_det - psi_coef_tmp(i,1) = psi_coef(i,igood_state_1) - psi_coef_tmp(i,2) = psi_coef(i,igood_state_2) - enddo - call save_wavefunction_general(N_det,2,psi_det,n_det,psi_coef_tmp) - deallocate(psi_coef_tmp) - - - -end