From cc71d6f83dcacf040863d9ba739887f7501af09d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 17 Feb 2016 12:22:45 +0100 Subject: [PATCH] Added OVB_effective_Hamiltonian --- .../Dressed_Ref_Hamiltonian.main.irp.f | 37 ++++++ .../NEEDED_CHILDREN_MODULES | 1 + plugins/Dressed_Ref_Hamiltonian/README.rst | 16 +++ .../dressed_eigenvectors.irp.f | 41 +++++++ .../dressed_hamiltonian.irp.f | 46 ++++++++ .../print_CAS_effective_Hamiltonian.irp.f | 108 ++++++++++++++++++ plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES | 2 +- plugins/MRCC_Utils/mrcc_dummy.irp.f | 4 + .../NEEDED_CHILDREN_MODULES | 1 + .../OVB_effective_H.irp.f | 59 ++++++++++ plugins/OVB_effective_Hamiltonian/README.rst | 13 +++ .../print_OVB_effective_H_diagonalized.irp.f | 101 ++++++++++++++++ ...save_wf_only_ionic_and_1p_amplitudes.irp.f | 90 +++++++++++++++ ...ve_wf_only_neutral_and_1p_amplitudes.irp.f | 88 ++++++++++++++ 14 files changed, 606 insertions(+), 1 deletion(-) create mode 100644 plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f create mode 100644 plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Dressed_Ref_Hamiltonian/README.rst create mode 100644 plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f create mode 100644 plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f create mode 100644 plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f create mode 100644 plugins/MRCC_Utils/mrcc_dummy.irp.f create mode 100644 plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES create mode 100644 plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f create mode 100644 plugins/OVB_effective_Hamiltonian/README.rst create mode 100644 plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f create mode 100644 plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f create mode 100644 plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f diff --git a/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f b/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f new file mode 100644 index 00000000..2e25f431 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f @@ -0,0 +1,37 @@ +program Dressed_Ref_Hamiltonian implicit none + BEGIN_DOC +! TODO + END_DOC + print *, ' _/ ' + print *, ' -:\_?, _Jm####La ' + print *, 'J"(:" > _]#AZ#Z#UUZ##, ' + print *, '_,::./ %(|i%12XmX1*1XL _?, ' + print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( ' + print *, ' .:< ]J=mQD?WXn|,)nr" ' + print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" ' + print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 ' + print *, ' "23Z#1S2oo2XXSnnnoSo2>v" ' + print *, ' miX#L -~`""!!1}oSoe|i7 ' + print *, ' 4cn#m, v221=|v[ ' + print *, ' ]hI3Zma,;..__wXSe=+vo ' + print *, ' ]Zov*XSUXXZXZXSe||vo2 ' + print *, ' ]Z#>=|< ' + print *, ' -ziiiii||||||+||==+> ' + print *, ' -%|+++||=|=+|=|==/ ' + print *, ' -a>====+|====-:- ' + print *, ' "~,- -- /- ' + print *, ' -. )> ' + print *, ' .~ +- ' + print *, ' . .... : . ' + print *, ' -------~ ' + print *, '' +end diff --git a/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES b/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..55c429bc --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +MRCC_Utils diff --git a/plugins/Dressed_Ref_Hamiltonian/README.rst b/plugins/Dressed_Ref_Hamiltonian/README.rst new file mode 100644 index 00000000..71f2f099 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/README.rst @@ -0,0 +1,16 @@ +======================= +Dressed_Ref_Hamiltonian +======================= +The following modules proposes to build an effective Hamiltonian +spanned on the reference determinants supposed to be the CAS ones. +The effective matrix Hamiltonian are built using the multi parentage +proposal used in the MR-CCSD formalism of Giner et. al. (JCP, 144, 064101 (2016); doi: 10.1063/1.4940781) + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f b/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f new file mode 100644 index 00000000..2e66ff16 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f @@ -0,0 +1,41 @@ + BEGIN_PROVIDER [double precision, psi_ref_coef_dressed, (n_det_ref,N_states) ] +&BEGIN_PROVIDER [double precision, energies_ref_dressed, (N_states) ] + implicit none + integer :: i,j,k,l,istate,igoodstate + double precision, allocatable :: H_matrix_tmp(:,:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:),psi_coef_ref_tmp(:) + double precision :: accu, accu1 + allocate(H_matrix_tmp(n_det_ref,n_det_ref)) + allocate(eigvalues(n_det_ref)) + allocate(eigvectors(n_det_ref,n_det_ref)) + allocate(psi_coef_ref_tmp(n_det_ref)) + do istate = 1, N_states + accu1 = 0.d0 + do j = 1, n_det_ref + accu1 += psi_ref_coef(j,istate)**2 ! norm of the "istate" eigenvector in the projected in the reference space + do k = 1, n_det_ref + H_matrix_tmp(j,k) = hamiltonian_total_dressed(j,k,istate) + enddo + enddo + accu1 = 1.d0/dsqrt(accu1) + do j = 1, n_det_ref + psi_coef_ref_tmp(j) = psi_ref_coef(j,istate) * accu1 + enddo + call lapack_diagd(eigvalues,eigvectors,H_matrix_tmp,n_det_ref,n_det_ref) + do j = 1, n_det_ref + accu = 0.d0 + do k = 1, n_det_ref + accu += eigvectors(k,j) * psi_coef_ref_tmp(k) + enddo + if(dabs(accu).gt.0.9d0)then + igoodstate = j + exit + endif + enddo + energies_ref_dressed(istate) = eigvalues(igoodstate) + do j = 1,n_det_ref + psi_ref_coef_dressed(j,istate) = eigvectors(j,igoodstate) + enddo + enddo + +END_PROVIDER diff --git a/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f b/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f new file mode 100644 index 00000000..90b27e07 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f @@ -0,0 +1,46 @@ +BEGIN_PROVIDER [double precision, dressing_ref_hamiltonian, (n_det_ref,n_det_ref,N_states)] + implicit none + integer :: i,j,k,l + integer :: ii,jj,istate + double precision :: hij,sec_order,H_ref(N_det_ref),hik,hkl + integer :: idx(0:N_det_ref) + double precision :: accu_negative,accu_positive,phase + integer :: degree_exc_ionic,degree_exc_neutral,exc(0:2,2,2) + dressing_ref_hamiltonian = 0.d0 + accu_negative = 0.d0 + accu_positive = 0.d0 + integer :: h1,p1,h2,p2,s1,s2 + do istate = 1, N_states + do i = 1, N_det_non_ref + call filter_connected_i_H_psi0(psi_ref,psi_non_ref(1,1,i),N_int,N_det_ref,idx) + H_ref = 0.d0 + do ii=1,idx(0) + k = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(psi_ref(1,1,k),psi_non_ref(1,1,i),N_int,hij) + H_ref(k) = hij + enddo + do ii= 1, idx(0) + k = idx(ii) + hik = H_ref(k) * lambda_mrcc(istate,i) + do jj = 1, idx(0) + l = idx(jj) + dressing_ref_hamiltonian(k,l,istate) += hik * H_ref(l) + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, hamiltonian_total_dressed, (n_det_ref,n_det_ref,N_states)] + implicit none + integer :: i,j,k + do k = 1, N_states + do i = 1, N_det_ref + do j = 1, N_det_ref + hamiltonian_total_dressed(j,i,k) = dressing_ref_hamiltonian(j,i,k) + ref_hamiltonian_matrix(j,i) + enddo + enddo + enddo + +END_PROVIDER diff --git a/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f b/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f new file mode 100644 index 00000000..aab716ac --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f @@ -0,0 +1,108 @@ +program print + read_wf = .True. + touch read_wf + call provide_all_stuffs +end +subroutine provide_all_stuffs + implicit none + provide ref_hamiltonian_matrix dressing_ref_hamiltonian + integer :: i,j,istate + double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + double precision, allocatable :: H_dressed(:,:) + double precision, allocatable :: H_print(:,:) + double precision :: accu_norm + allocate (H_dressed(N_det_ref,N_det_ref)) + allocate (H_print(N_det_ref,N_det_ref)) + allocate (psi_restart_ref_normalized(N_det_ref)) + allocate (psi_ref_zeroth_order(N_det_ref)) + print*,'# nuclear_repulsion = ',nuclear_repulsion + allocate (psi_ref_dressed(N_det_ref)) + allocate (eigvalues(N_det_ref)) + allocate (eigvectors(N_det_ref,N_det_ref)) + + + + do istate= 1, N_states + do i = 1, N_det_ref + do j = 1, N_det_ref + H_print(i,j) = ref_hamiltonian_matrix(j,i) + enddo + enddo + do i = 1, N_det_ref + H_print(i,i) -= ref_hamiltonian_matrix(1,1) + enddo + print*,'Ref Hamiltonian matrix emelent = ',ref_hamiltonian_matrix(1,1) + print*,'ISTATE = ',istate + accu_norm = 0.d0 + do i = 1, N_det_ref + accu_norm += psi_ref_coef(i,1) * psi_ref_coef(i,1) + enddo + print*,'accu_norm = ',accu_norm + accu_norm = 1.d0/dsqrt(accu_norm) + do i = 1, N_det_ref + psi_restart_ref_normalized(i) = psi_ref_coef(i,istate)* accu_norm + enddo + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') H_print(i,:) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX DRESSING' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') dressing_ref_hamiltonian(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + do i = 1, N_det_ref + do j = 1, N_det_ref + H_dressed(j,i) = ref_hamiltonian_matrix(j,i) + dressing_ref_hamiltonian(j,i,istate) + H_print(i,j) += dressing_ref_hamiltonian(j,i,istate) + enddo + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'TOTAL DRESSED H MATRIX ' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') H_print(i,:) + enddo + print*,'' + print*,'' + print*,'' + + + call lapack_diagd(eigvalues,eigvectors,ref_hamiltonian_matrix,n_det_ref,n_det_ref) + do i = 1, N_det_ref + psi_ref_zeroth_order(i) = eigvectors(i,istate) + enddo + + + call lapack_diagd(eigvalues,eigvectors,H_dressed,n_det_ref,n_det_ref) + do i = 1, N_det_ref + psi_ref_dressed(i) = eigvectors(i,istate) + enddo + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = 1, N_det_ref + write(*,'(10(F10.7 ,4X))') psi_ref_coef(i,istate)/psi_ref_coef(1,istate), psi_ref_dressed(i)/psi_ref_dressed(1),psi_ref_zeroth_order(i)/psi_ref_zeroth_order(1) + enddo + enddo + + deallocate (H_dressed) + deallocate (H_print) + deallocate (psi_restart_ref_normalized) + deallocate (psi_ref_zeroth_order) + deallocate (psi_ref_dressed) + + deallocate (eigvalues) + deallocate (eigvectors) + +end diff --git a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES index 5b16423e..7392852a 100644 --- a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES +++ b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full Psiref_Utils +Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS diff --git a/plugins/MRCC_Utils/mrcc_dummy.irp.f b/plugins/MRCC_Utils/mrcc_dummy.irp.f new file mode 100644 index 00000000..8f1deda8 --- /dev/null +++ b/plugins/MRCC_Utils/mrcc_dummy.irp.f @@ -0,0 +1,4 @@ +program pouet + + +end diff --git a/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES b/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..01b1f4b9 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Dressed_Ref_Hamiltonian OVB diff --git a/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f b/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f new file mode 100644 index 00000000..4ef79aeb --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f @@ -0,0 +1,59 @@ + BEGIN_PROVIDER [double precision, H_OVB_dressing, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of all the + END_DOC + implicit none + integer :: i,j,istate,k,l + double precision :: accu,hij + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + accu = 0.d0 + do istate = 1, N_states + do k = 1, ionic_index(i,0) + do l = 1, ionic_index(j,0) + hij = dressing_ref_hamiltonian(ionic_index(i,k),ionic_index(j,l),istate) +! accu += psi_ref_coef(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * & +! psi_ref_coef(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + accu += psi_ref_coef_dressed(ionic_index(i,k),istate) * normalization_factor_ionic_dressed(i,istate) * & + psi_ref_coef_dressed(ionic_index(j,l),istate) * normalization_factor_ionic_dressed(j,istate) * hij + enddo + enddo + H_OVB_dressing(i,j,istate) = accu + enddo + enddo + enddo + END_PROVIDER + + + + BEGIN_PROVIDER [double precision, H_OVB_total_dressed, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of all the + END_DOC + implicit none + integer :: i,j,istate + double precision :: accu,hij + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + do istate = 1, N_states + H_OVB_total_dressed(i,j,istate) = H_OVB_dressing(i,j,istate) + H_OVB_naked(i,j,istate) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, normalization_factor_ionic_dressed, (min_number_ionic:max_number_ionic, N_states) ] + implicit none + integer :: i,j,istate,k + double precision :: accu + do j = min_number_ionic, max_number_ionic + do istate = 1, N_states + accu = 0.d0 + do k = 1, ionic_index(j,0) + accu += psi_ref_coef_dressed(ionic_index(j,k),istate) **2 + enddo + normalization_factor_ionic_dressed(j,istate) = 1.d0/dsqrt(accu) + enddo + enddo + + END_PROVIDER diff --git a/plugins/OVB_effective_Hamiltonian/README.rst b/plugins/OVB_effective_Hamiltonian/README.rst new file mode 100644 index 00000000..a21ffcc6 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/README.rst @@ -0,0 +1,13 @@ +========================= +OVB_effective_Hamiltonian +========================= +Dressing of the OVB matrix by use of the Dressed_Ref_Hamiltonian dressing + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f b/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f new file mode 100644 index 00000000..bbd29b8e --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f @@ -0,0 +1,101 @@ +program print + read_wf = .True. + touch read_wf + call provide_all_stuffs +end +subroutine provide_all_stuffs + implicit none + provide ref_hamiltonian_matrix dressing_ref_hamiltonian + integer :: i,j,istate + double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + double precision, allocatable :: H_naked(:,:) + double precision, allocatable :: H_dressed(:,:) + double precision, allocatable :: H_print(:,:) + double precision :: accu_norm + allocate (H_dressed(max_number_ionic+1,max_number_ionic+1)) + allocate (H_print(min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)) + allocate (H_naked(max_number_ionic+1,max_number_ionic+1)) + allocate (psi_restart_ref_normalized(min_number_ionic:max_number_ionic)) + allocate (psi_ref_zeroth_order(min_number_ionic:max_number_ionic)) + print*,'# nuclear_repulsion = ',nuclear_repulsion + allocate (psi_ref_dressed(min_number_ionic:max_number_ionic)) + allocate (eigvalues(max_number_ionic+1)) + allocate (eigvectors(max_number_ionic+1,max_number_ionic+1)) + + do istate= 1, N_states + print*,'ISTATE = ',istate + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_print(i,j) = H_OVB_naked(j,i,istate) + enddo + enddo + do i = min_number_ionic,max_number_ionic + H_print(i,i) -= H_OVB_naked(min_number_ionic,min_number_ionic,istate) + enddo + + print*,'Ref Hamiltonian matrix emelent = ',H_OVB_naked(min_number_ionic,min_number_ionic,istate) + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:) + enddo + print*,'CAS MATRIX DRESSING' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_dressing(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX DRESSED ' + print*,'' + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_print(i,j) = H_OVB_total_dressed(j,i,istate) + enddo + enddo + do i = min_number_ionic,max_number_ionic + H_print(i,i) -= H_OVB_total_dressed(min_number_ionic,min_number_ionic,istate) + enddo + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:) + enddo + print*,'' + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_dressed(j+1,i+1) = H_OVB_total_dressed(i,j,istate) + H_naked(j+1,i+1) = H_OVB_naked(i,j,istate) + enddo + enddo + + call lapack_diagd(eigvalues,eigvectors,H_naked,max_number_ionic+1,max_number_ionic+1) + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = min_number_ionic,max_number_ionic + psi_ref_zeroth_order(i) = eigvectors(i+1,istate) + enddo + + + + call lapack_diagd(eigvalues,eigvectors,H_dressed,max_number_ionic+1,max_number_ionic+1) + do i = min_number_ionic,max_number_ionic + psi_ref_dressed(i) = eigvectors(i+1,istate) + enddo + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = min_number_ionic,max_number_ionic + write(*,'(10(F10.7 ,4X))') psi_ref_dressed(i)/psi_ref_dressed(min_number_ionic) ,psi_ref_zeroth_order(i)/psi_ref_zeroth_order(min_number_ionic) + enddo + enddo + + deallocate (H_dressed) + deallocate (H_naked) + deallocate (psi_restart_ref_normalized) + deallocate (psi_ref_zeroth_order) + deallocate (psi_ref_dressed) + + deallocate (eigvalues) + deallocate (eigvectors) + +end diff --git a/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f b/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f new file mode 100644 index 00000000..4398b152 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f @@ -0,0 +1,90 @@ +program save_wf + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + use bitmasks + integer :: i,j,k,l + integer(bit_kind), allocatable :: psi_save_final(:,:,:) + double precision, allocatable :: psi_coef_save_final(:,:) + integer :: index_ref_determinants_save(psi_det_size) + integer :: n_det_ref_determinants_save + integer :: index_non_ref_determinants_save(psi_det_size) + integer :: n_det_non_ref_determinants_save + + integer :: n_det_save_final + integer :: number_of_particles + n_det_ref_determinants_save = 0 + integer :: ionic_level + ionic_level = 1 + do i = 1, ionic_index(ionic_level,0) ! number of determinants in the ref wf that are neutrals + n_det_ref_determinants_save +=1 + index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(ionic_level,i) + enddo + ! save all the 1p determinants in order to have the single excitations + ! on the top of the neutral structures + n_det_non_ref_determinants_save = 0 + do i = 1, N_det_non_ref + if(number_of_particles(psi_non_ref(1,1,i))==1)then + n_det_non_ref_determinants_save +=1 + index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i + endif + enddo + print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save + print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save + n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save + allocate (psi_save_final(N_int,2,n_det_save_final)) + allocate (psi_coef_save_final(n_det_save_final,1)) + integer :: n_det_tmp + n_det_tmp = 0 + do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i)) + enddo + psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1) + enddo + pause + do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i)) + enddo + accu = 0.d0 + double precision :: t_ik,hij + do j = 1, n_det_ref_determinants_save + call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij) + t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i)) + accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik + enddo + psi_coef_save_final(n_det_tmp,1) = accu + enddo + double precision :: accu + accu = 0.d0 + do i = 1, n_det_save_final + accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1) + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, n_det_save_final + psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1) + enddo + + do i = 1, n_det_save_final + print*,'' + print*,'Det' + call debug_det(psi_save_final(1,1,i),N_int) + print*,'coef = ',psi_coef_save_final(i,1) + enddo + + call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final) + deallocate (psi_save_final) + deallocate (psi_coef_save_final) + + +end diff --git a/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f b/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f new file mode 100644 index 00000000..a4b83d33 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f @@ -0,0 +1,88 @@ +program save_wf + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + use bitmasks + integer :: i,j,k,l + integer(bit_kind), allocatable :: psi_save_final(:,:,:) + double precision, allocatable :: psi_coef_save_final(:,:) + integer :: index_ref_determinants_save(psi_det_size) + integer :: n_det_ref_determinants_save + integer :: index_non_ref_determinants_save(psi_det_size) + integer :: n_det_non_ref_determinants_save + + integer :: n_det_save_final + integer :: number_of_particles + n_det_ref_determinants_save = 0 + do i = 1, ionic_index(0,0) ! number of determinants in the ref wf that are neutrals + n_det_ref_determinants_save +=1 + index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(0,i) + enddo + ! save all the 1p determinants in order to have the single excitations + ! on the top of the neutral structures + n_det_non_ref_determinants_save = 0 + do i = 1, N_det_non_ref + if(number_of_particles(psi_non_ref(1,1,i))==1)then + n_det_non_ref_determinants_save +=1 + index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i + endif + enddo + print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save + print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save + n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save + allocate (psi_save_final(N_int,2,n_det_save_final)) + allocate (psi_coef_save_final(n_det_save_final,1)) + integer :: n_det_tmp + n_det_tmp = 0 + do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i)) + enddo + psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1) + enddo + pause + do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i)) + enddo + accu = 0.d0 + double precision :: t_ik,hij + do j = 1, n_det_ref_determinants_save + call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij) + t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i)) + accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik + enddo + psi_coef_save_final(n_det_tmp,1) = accu + enddo + double precision :: accu + accu = 0.d0 + do i = 1, n_det_save_final + accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1) + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, n_det_save_final + psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1) + enddo + + do i = 1, n_det_save_final + print*,'' + print*,'Det' + call debug_det(psi_save_final(1,1,i),N_int) + print*,'coef = ',psi_coef_save_final(i,1) + enddo + + call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final) + deallocate (psi_save_final) + deallocate (psi_coef_save_final) + + +end