From 27a121bc1ba960bf1e778af2388e8c46b756e9e4 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 17 Feb 2016 10:52:57 +0100 Subject: [PATCH] added print_hcc and OVB plugin --- config/ifort.cfg | 2 +- plugins/Hartree_Fock/EZFIO.cfg | 7 + plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES | 2 +- plugins/Hartree_Fock/diagonalize_fock.irp.f | 74 +-- plugins/OVB/NEEDED_CHILDREN_MODULES | 1 + plugins/OVB/README.rst | 20 + plugins/OVB/ovb_components.irp.f | 510 +++++++++++++++++++ plugins/OVB/print_ovb.irp.f | 27 + plugins/Properties/print_hcc.irp.f | 17 + plugins/Psiref_Utils/psi_ref_utils.irp.f | 21 +- src/Bitmask/bitmasks.irp.f | 61 ++- src/Determinants/usefull_for_ovb.irp.f | 283 ++++++++++ 12 files changed, 967 insertions(+), 58 deletions(-) create mode 100644 plugins/OVB/NEEDED_CHILDREN_MODULES create mode 100644 plugins/OVB/README.rst create mode 100644 plugins/OVB/ovb_components.irp.f create mode 100644 plugins/OVB/print_ovb.irp.f create mode 100644 plugins/Properties/print_hcc.irp.f create mode 100644 src/Determinants/usefull_for_ovb.irp.f diff --git a/config/ifort.cfg b/config/ifort.cfg index f8c36e4e..cc848cba 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index d8207cc4..2fa29cf0 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -26,3 +26,10 @@ default: Huckel type: double precision doc: Calculated HF energy interface: ezfio + +[no_oa_or_av_opt] +type: logical +doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the SCF procedure +interface: ezfio,provider,ocaml +default: False + diff --git a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES index 85bdd3ad..6fb87e35 100644 --- a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES +++ b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Integrals_Bielec MOGuess +Integrals_Bielec MOGuess Bitmask diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f index 7c424aeb..c80077b3 100644 --- a/plugins/Hartree_Fock/diagonalize_fock.irp.f +++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f @@ -11,63 +11,35 @@ double precision, allocatable :: work(:), F(:,:), S(:,:) -! if (mo_tot_num == ao_num) then -! ! Solve H.C = E.S.C in AO basis set -! -! allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) -! do j=1,ao_num -! do i=1,ao_num -! S(i,j) = ao_overlap(i,j) -! F(i,j) = Fock_matrix_ao(i,j) -! enddo -! enddo -! -! n = ao_num -! lwork = 1+6*n + 2*n*n -! liwork = 3 + 5*n -! -! allocate(work(lwork), iwork(liwork) ) -! -! lwork = -1 -! liwork = -1 -! -! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! -! if (info /= 0) then -! print *, irp_here//' failed : ', info -! stop 1 -! endif -! lwork = int(work(1)) -! liwork = iwork(1) -! deallocate(work,iwork) -! allocate(work(lwork), iwork(liwork) ) -! -! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! -! if (info /= 0) then -! print *, irp_here//' failed : ', info -! stop 1 -! endif -! do j=1,mo_tot_num -! do i=1,ao_num -! eigenvectors_Fock_matrix_mo(i,j) = F(i,j) -! enddo -! enddo -! -! deallocate(work, iwork, F, S) -! -! else -! - ! Solve H.C = E.C in MO basis set - allocate( F(mo_tot_num_align,mo_tot_num) ) do j=1,mo_tot_num do i=1,mo_tot_num F(i,j) = Fock_matrix_mo(i,j) enddo enddo + if(no_oa_or_av_opt)then + integer :: iorb,jorb + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + enddo + endif + + ! Insert level shift here diff --git a/plugins/OVB/NEEDED_CHILDREN_MODULES b/plugins/OVB/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..05706cac --- /dev/null +++ b/plugins/OVB/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants Psiref_CAS diff --git a/plugins/OVB/README.rst b/plugins/OVB/README.rst new file mode 100644 index 00000000..13488c7e --- /dev/null +++ b/plugins/OVB/README.rst @@ -0,0 +1,20 @@ +======================= +OVB +======================= +The present module proposes an orthogonal Valence Bond analysis +of the wave function, that are the printing of the various Hamiltonian +matrix elements on the basis of the level of ionicity of the components +of the wave function. + +Assumptions : it supposes that you have some orthogonal local orbitals within +the active space and that you performed a CI within the active orbitals. +Such CI might be complete or not, no matter. + +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/ovb_components.irp.f b/plugins/OVB/ovb_components.irp.f new file mode 100644 index 00000000..82f0f931 --- /dev/null +++ b/plugins/OVB/ovb_components.irp.f @@ -0,0 +1,510 @@ + + use bitmasks + BEGIN_PROVIDER [integer, max_number_ionic] +&BEGIN_PROVIDER [integer, min_number_ionic] + BEGIN_DOC + ! Maximum and minimum number of ionization in psi_ref + END_DOC + implicit none + integer :: i,j + integer :: n_closed_shell_cas + max_number_ionic = 0 + min_number_ionic = 100000 + do i = 1, N_det_ref + j = n_closed_shell_cas(psi_ref(1,1,i),n_int) + if(j> max_number_ionic)then + max_number_ionic = j + endif + if(j< min_number_ionic)then + min_number_ionic = j + endif + + enddo + print*,'max_number_ionic = ',max_number_ionic + print*,'min_number_ionic = ',min_number_ionic +END_PROVIDER + + BEGIN_PROVIDER [integer, ionic_index, (min_number_ionic:max_number_ionic,0:N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_ionic, (min_number_ionic:max_number_ionic, N_states)] + BEGIN_DOC + ! Index of the various determinants in psi_ref according to their level of ionicity + ! ionic_index(i,0) = number of determinants in psi_ref having the degree of ionicity "i" + ! ionic_index(i,j) = index of the determinants having the degree of ionicity "i" + END_DOC + implicit none + integer :: i,j,k + integer :: n_closed_shell_cas + double precision :: accu + ionic_index = 0 + do i = 1, N_det_ref + j = n_closed_shell_cas(psi_ref(1,1,i),n_int) + ionic_index(j,0) +=1 + ionic_index(j,ionic_index(j,0)) = i + enddo + do i = min_number_ionic,max_number_ionic + accu = 0.d0 + do j = 1, N_states + do k = 1, ionic_index(i,0) + accu += psi_ref_coef_diagonalized(ionic_index(i,k),j) * psi_ref_coef_diagonalized(ionic_index(i,k),j) + enddo + normalization_factor_ionic(i,j) = 1.d0/dsqrt(accu) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [double precision, H_OVB_naked, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of contracted forms in terms of ionic structures + 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 + do istate = 1, N_states + accu = 0.d0 + do k = 1, ionic_index(i,0) + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + enddo + H_OVB_naked(i,j,istate) = accu + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [integer, n_couples_act_orb] + implicit none + n_couples_act_orb = 3 + END_PROVIDER + + BEGIN_PROVIDER [integer, couples_act_orb, (n_couples_act_orb,2) ] + implicit none + + couples_act_orb(1,1) = 20 + couples_act_orb(1,2) = 21 + couples_act_orb(2,1) = 22 + couples_act_orb(2,2) = 23 + couples_act_orb(3,1) = 24 + couples_act_orb(3,2) = 25 + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom , (n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC +! Hamiltonian matrix elements between the various contracted functions +! that have a negative charge on a given active orbital + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + do i = 1, n_act_orb + j = i ! Diagonal part + norm = 0.d0 + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + norm += psi_coef_mono_ionic_on_given_atom(k,i) **2 + do ii = 1, N_int + key_1(ii,1) = psi_det_mono_ionic_on_given_atom(ii,1,k,i) + key_1(ii,2) = psi_det_mono_ionic_on_given_atom(ii,2,k,i) + enddo + do l = 1, n_det_ionic_on_given_atom(j) + do jj = 1, N_int + key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j) + key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom(i,j) = accu + + + do j = i+1, n_act_orb ! Extra diagonal part + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + do jj = 1, N_int + key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i) + key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i) + enddo + do l = 1, n_det_ionic_on_given_atom(j) + do jj = 1, N_int + key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j) + key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom(i,j) = accu + H_matrix_between_ionic_on_given_atom(j,i) = accu + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom_and_others , (n_act_orb,min_number_ionic:max_number_ionic)] + implicit none + use bitmasks + BEGIN_DOC +! Hamiltonian matrix elements between the various contracted functions +! that have a negative charge on a given active orbital +! and all the other fully contracted OVB structures + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + + do i = 1, n_act_orb + do j = min_number_ionic,max_number_ionic + if(j==1)then + H_matrix_between_ionic_on_given_atom_and_others(i,j) = 0.d0 + endif + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + do jj = 1, N_int + key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i) + key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i) + enddo + do l = 1, ionic_index(j,0) + do ii = 1, N_int + key_2(ii,1) = psi_det_ovb(ii,1,l,j) + key_2(ii,2) = psi_det_ovb(ii,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_ovb(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom_and_others(i,j) = accu + enddo + enddo + + print*,'H_matrix_between_ionic_on_given_atom_and_others' + print*,'' + do i = 1, n_act_orb + write(*,'(I3,X,100(F16.7))'),H_matrix_between_ionic_on_given_atom_and_others(i,:) + enddo + + + +END_PROVIDER + + BEGIN_PROVIDER [integer, n_det_ionic_on_given_atom, (n_act_orb)] +&BEGIN_PROVIDER [double precision, normalization_factor_ionic_on_given_atom, (n_act_orb) ] +&BEGIN_PROVIDER [double precision, psi_coef_mono_ionic_on_given_atom, (N_det_ref,n_act_orb) ] +&BEGIN_PROVIDER [integer(bit_kind), psi_det_mono_ionic_on_given_atom, (N_int,2,N_det_ref,n_act_orb)] + implicit none + use bitmasks + BEGIN_DOC +! number of determinants that are mono ionic with the negative charge +! on a given atom, normalization_factor, array of determinants,and coefficients + END_DOC + integer :: i,j,k,l + ionicity_level = 1 + integer :: ionicity_level + logical :: doubly_occupied_array(n_act_orb) + n_det_ionic_on_given_atom = 0 + normalization_factor_ionic_on_given_atom = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + call give_index_of_doubly_occ_in_active_space(psi_det(1,1,ionic_index(ionicity_level,i)),doubly_occupied_array) + do j = 1, n_act_orb + if(doubly_occupied_array(j))then + n_det_ionic_on_given_atom(j) += 1 + normalization_factor_ionic_on_given_atom(j) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + do k = 1, N_int + psi_det_mono_ionic_on_given_atom(k,1,n_det_ionic_on_given_atom(j),j) = psi_det(k,1,ionic_index(ionicity_level,i)) + psi_det_mono_ionic_on_given_atom(k,2,n_det_ionic_on_given_atom(j),j) = psi_det(k,2,ionic_index(ionicity_level,i)) + enddo + psi_coef_mono_ionic_on_given_atom(n_det_ionic_on_given_atom(j),j) = psi_ref_coef_diagonalized(ionic_index(1,i),1) + endif + enddo + enddo + integer :: i_count + i_count = 0 + do j = 1, n_act_orb + i_count += n_det_ionic_on_given_atom(j) + normalization_factor_ionic_on_given_atom(j) = 1.d0/dsqrt(normalization_factor_ionic_on_given_atom(j)) + enddo + if(i_count.ne.ionic_index(ionicity_level,0))then + print*,'PB with n_det_ionic_on_given_atom' + print*,'i_count = ',i_count + print*,'ionic_index(ionicity_level,0)',ionic_index(ionicity_level,0) + stop + endif + do j = 1, n_act_orb + do i = 1, n_det_ionic_on_given_atom(j) + psi_coef_mono_ionic_on_given_atom(i,j) = psi_coef_mono_ionic_on_given_atom(i,j) * normalization_factor_ionic_on_given_atom(j) + enddo + enddo + + + END_PROVIDER + + BEGIN_PROVIDER [integer(bit_kind), psi_det_ovb, (N_int,2,N_det_ref,min_number_ionic:max_number_ionic)] +&BEGIN_PROVIDER [double precision, psi_coef_ovb, (N_det_ref,min_number_ionic:max_number_ionic) ] + implicit none + BEGIN_DOC +! Array of the determinants belonging to each ovb structures (neutral, mono ionic, bi ionic etc ...) +! together with the arrays of coefficients + END_DOC + integer :: i,j,k,l + use bitmasks + integer :: ionicity_level,i_count + double precision :: accu + + do ionicity_level = min_number_ionic,max_number_ionic + accu = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + psi_det_ovb(j,1,i,ionicity_level) = psi_det(j,1,ionic_index(ionicity_level,i)) + psi_det_ovb(j,2,i,ionicity_level) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + psi_coef_ovb(i,ionicity_level) = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) * normalization_factor_ionic(ionicity_level,1) + accu += psi_coef_ovb(i,ionicity_level)**2 + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, ionic_index(ionicity_level,0) + psi_coef_ovb(i,ionicity_level) = psi_coef_ovb(i,ionicity_level) * accu + enddo + accu = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + accu += psi_coef_ovb(i,ionicity_level) **2 + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, H_matrix_psi_det_ovb, (min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)] + implicit none + BEGIN_DOC +! H matrix between the fully contracted OVB forms + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + use bitmasks + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + accu = 0.d0 + do k = 1, ionic_index(i,0) + do ii = 1, N_int + key_1(ii,1) = psi_det_ovb(ii,1,k,i) + key_1(ii,2) = psi_det_ovb(ii,2,k,i) + enddo + do l = 1, ionic_index(j,0) + do ii = 1, N_int + key_2(ii,1) = psi_det_ovb(ii,1,l,j) + key_2(ii,2) = psi_det_ovb(ii,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_ovb(l,j) * psi_coef_ovb(k,i) * hij + enddo + enddo + H_matrix_psi_det_ovb(i,j) = accu + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [integer, number_first_ionic_couples] +&BEGIN_PROVIDER [logical , is_a_first_ionic_couple, (N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_special_first_ionic, (2)] + implicit none + BEGIN_DOC + ! Number of determinants belonging to the class of first ionic + ! AND that have a couple of positive/negative charge belonging + ! to a couple of orbital couples_act_orb + ! If is_a_first_ionic_couple(i) = .True. then this determinant is a first ionic + ! and have a couple of positive/negative charge belonging + ! to a couple of orbital couples_act_orb + ! normalization factor (1) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .True.) + ! normalization factor (2) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .False.) + END_DOC + integer :: i,j + use bitmasks + number_first_ionic_couples = 0 + integer :: ionicity_level + logical :: couples_out(0:n_couples_act_orb) + integer(bit_kind) :: key_tmp(N_int,2) + ionicity_level = 1 + normalization_factor_special_first_ionic = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i)) + key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + call doubly_occ_empty_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out) + if(couples_out(0))then + number_first_ionic_couples +=1 + is_a_first_ionic_couple(i) = .True. + normalization_factor_special_first_ionic(1) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + else + is_a_first_ionic_couple(i) = .False. + normalization_factor_special_first_ionic(2) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + endif + enddo + normalization_factor_special_first_ionic(1) = 1.d0/dsqrt(normalization_factor_special_first_ionic(1)) + normalization_factor_special_first_ionic(2) = 1.d0/dsqrt(normalization_factor_special_first_ionic(2)) + print*,'number_first_ionic_couples = ',number_first_ionic_couples + END_PROVIDER + + + BEGIN_PROVIDER [integer, number_neutral_no_hund_couples] +&BEGIN_PROVIDER [logical , is_a_neutral_no_hund_couple, (N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_neutra_no_hund_couple, (2)] +&BEGIN_PROVIDER [double precision, ratio_hund_no_hund ] + implicit none + BEGIN_DOC + ! Number of determinants belonging to the class of neutral determinants + ! AND that have a couple of alpha beta electrons in couple of orbital couples_act_orb + ! If is_a_neutral_no_hund_couple(i) = .True. then this determinant is a neutral determinants + ! and have a a couple of alpha beta electrons in couple of orbital couples_act_orb + ! normalization factor (1) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .True.) + ! normalization factor (2) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .False.) + END_DOC + integer :: i,j + use bitmasks + number_neutral_no_hund_couples = 0 + integer :: ionicity_level + logical :: couples_out(0:n_couples_act_orb) + integer(bit_kind) :: key_tmp(N_int,2) + integer :: ifirst_hund,ifirst_no_hund + double precision :: coef_ref_hund,coef_ref_no_hund + ifirst_hund = 0 + ifirst_no_hund = 0 + ionicity_level = 0 + normalization_factor_neutra_no_hund_couple = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i)) + key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + call neutral_no_hund_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out) + if(couples_out(0))then + if(ifirst_no_hund == 0)then + coef_ref_no_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) + ifirst_no_hund = 1 + endif + number_neutral_no_hund_couples +=1 + is_a_neutral_no_hund_couple(i) = .True. + normalization_factor_neutra_no_hund_couple(1) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2 + else + if(ifirst_hund == 0)then + coef_ref_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) + ifirst_hund = 1 + endif + is_a_neutral_no_hund_couple(i) = .False. + normalization_factor_neutra_no_hund_couple(2) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2 + endif + enddo + ratio_hund_no_hund = coef_ref_no_hund/coef_ref_hund + + normalization_factor_neutra_no_hund_couple(1) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(1)) + normalization_factor_neutra_no_hund_couple(2) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(2)) + print*,'number_neutral_no_hund_couples = ',number_neutral_no_hund_couples + END_PROVIDER + + BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic, (2,min_number_ionic:max_number_ionic,n_states)] +&BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic_between_ionic, (2,2,n_states)] + BEGIN_DOC + ! H_OVB_naked_first_ionic(1,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the contracted ith ionic form + ! if i == 1 not defined + ! H_OVB_naked_first_ionic(2,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False + ! and the contracted ith ionic form + ! if i == 1 not defined + ! H_OVB_naked_first_ionic_between_ionic(1,1) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the first ionic determinants belonging to is_a_first_ionic_couple = True + ! H_OVB_naked_first_ionic_between_ionic(1,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the first ionic determinants belonging to is_a_first_ionic_couple = False + ! H_OVB_naked_first_ionic_between_ionic(2,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False + ! and the first ionic determinants belonging to is_a_first_ionic_couple = False + END_DOC + implicit none + integer :: i,j,istate,k,l + double precision :: accu_1,accu_2,hij + H_OVB_naked_first_ionic = 0.d0 + H_OVB_naked_first_ionic_between_ionic = 0.d0 + i = 1 + do j = min_number_ionic,max_number_ionic + if(j==1)cycle + do istate = 1, N_states + accu_1 = 0.d0 + accu_2 = 0.d0 + do k = 1, ionic_index(i,0) + if(is_a_first_ionic_couple(k))then + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu_1 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(1) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + else + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu_2 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(2) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + endif + enddo + H_OVB_naked_first_ionic(1,j,istate) = accu_1 + H_OVB_naked_first_ionic(2,j,istate) = accu_2 + enddo + enddo + + + do istate = 1, N_states + accu_1 = 0.d0 + accu_2 = 0.d0 + integer :: i_count + i_count = 0 + do k = 1, ionic_index(1,0) + do l = 1, ionic_index(1,0) + hij = ref_hamiltonian_matrix(ionic_index(1,k),ionic_index(1,l)) + accu_1 = hij * psi_ref_coef_diagonalized(ionic_index(1,k),istate) * psi_ref_coef_diagonalized(ionic_index(1,l),istate) + if(is_a_first_ionic_couple(k).and. is_a_first_ionic_couple(l))then + H_OVB_naked_first_ionic_between_ionic(1,1,istate) += accu_1 * normalization_factor_special_first_ionic(1) **2 + elseif(is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then + i_count += 1 + H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * & + normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2) +! elseif(is_a_first_ionic_couple(l).and. .not.is_a_first_ionic_couple(k))then +! i_count += 1 +! H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * & +! normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2) + elseif(.not.is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then + H_OVB_naked_first_ionic_between_ionic(2,2,istate) += accu_1 * normalization_factor_special_first_ionic(2) **2 + endif + enddo + enddo + enddo + print*,'i_count = ',i_count + print*,'number_first_ionic_couples**2 = ',ionic_index(1,0) * number_first_ionic_couples + + double precision :: convert_hartree_ev + convert_hartree_ev = 27.211399d0 + print*,'Special H matrix' + do i = 1,2 + write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic(i,:,1) + enddo + + print*,'Special H matrix bis' + do i = 1,2 + write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic_between_ionic(i,:,1) + enddo + + + END_PROVIDER + diff --git a/plugins/OVB/print_ovb.irp.f b/plugins/OVB/print_ovb.irp.f new file mode 100644 index 00000000..9e333c15 --- /dev/null +++ b/plugins/OVB/print_ovb.irp.f @@ -0,0 +1,27 @@ +program print_OVB + implicit none + read_wf = .True. + call provide_all + +end + +subroutine provide_all + implicit none + integer :: i,j,k,l,istate + do istate= 1, N_states + print*,'-------------------' + print*,'ISTATE = ',istate + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_naked(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + enddo + + +end + diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f new file mode 100644 index 00000000..f0091e1e --- /dev/null +++ b/plugins/Properties/print_hcc.irp.f @@ -0,0 +1,17 @@ +program print_hcc + implicit none + read_wf = .True. + touch read_wf + call test +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/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f index f9cf1303..fb45b13d 100644 --- a/plugins/Psiref_Utils/psi_ref_utils.irp.f +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -125,7 +125,7 @@ BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)] enddo END_PROVIDER - BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)] + BEGIN_PROVIDER [double precision, psi_ref_coef_diagonalized, (N_det_ref,N_states)] &BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)] implicit none integer :: i,j @@ -137,9 +137,11 @@ END_PROVIDER do i = 1, N_states psi_ref_energy_diagonalized(i) = eigenvalues(i) do j = 1, N_det_ref - psi_coef_ref_diagonalized(j,i) = eigenvectors(j,i) + psi_ref_coef_diagonalized(j,i) = eigenvectors(j,i) enddo enddo + deallocate (eigenvectors) + deallocate (eigenvalues) END_PROVIDER @@ -264,3 +266,18 @@ integer function get_index_in_psi_ref_sorted_bit(key,Nint) end +BEGIN_PROVIDER [double precision, ref_hamiltonian_matrix, (n_det_ref,n_det_ref)] + BEGIN_DOC + ! H matrix in the Reference space + END_DOC + implicit none + integer :: i,j + double precision :: hij + do i = 1, N_det_ref + do j = 1, N_det_ref + call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij) + ref_hamiltonian_matrix(i,j) = hij + enddo + enddo +END_PROVIDER + diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 70a84a42..50a95312 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -289,7 +289,12 @@ END_PROVIDER &BEGIN_PROVIDER [ integer, n_virt_orb ] implicit none BEGIN_DOC - ! Bitmasks for the inactive orbitals that are excited in post CAS method + ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited + ! in post CAS methods + ! n_inact_orb : Number of inactive orbitals + ! virt_bitmask : Bitmaks of vritual orbitals which are supposed to be recieve electrons + ! in post CAS methods + ! n_virt_orb : Number of virtual orbitals END_DOC logical :: exists integer :: j,i @@ -327,8 +332,14 @@ END_PROVIDER - BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] + BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] &BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)] + BEGIN_DOC + ! list_inact : List of the inactive orbitals which are supposed to be doubly excited + ! in post CAS methods + ! list_virt : List of vritual orbitals which are supposed to be recieve electrons + ! in post CAS methods + END_DOC implicit none integer :: occ_inact(N_int*bit_kind_size) integer :: itest,i @@ -348,6 +359,21 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1)) + reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2)) + enddo + END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] implicit none BEGIN_DOC @@ -376,7 +402,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)] implicit none BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks + ! Bitmask of the core orbitals that are never excited in post CAS method END_DOC integer :: i,j do i = 1, N_int @@ -385,6 +411,35 @@ END_PROVIDER enddo END_PROVIDER + BEGIN_PROVIDER [integer, list_core, (n_core_orb)] + BEGIN_DOC + ! List of the core orbitals that are never excited in post CAS method + END_DOC + implicit none + integer :: occ_core(N_int*bit_kind_size) + integer :: itest,i + occ_core = 0 + call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int) + ASSERT(itest==n_core_orb) + do i = 1, n_core_orb + list_core(i) = occ_core(i) + enddo + END_PROVIDER + + BEGIN_PROVIDER [ integer, n_core_orb ] + implicit none + BEGIN_DOC + ! Number of core orbitals that are never excited in post CAS method + END_DOC + logical :: exists + integer :: j,i + integer :: i_hole,i_part,i_gen + + n_core_orb = 0 + do j = 1, N_int + n_core_orb += popcnt(core_bitmask(j,1)) + enddo + END_PROVIDER BEGIN_PROVIDER [ integer, i_bitmask_gen ] diff --git a/src/Determinants/usefull_for_ovb.irp.f b/src/Determinants/usefull_for_ovb.irp.f new file mode 100644 index 00000000..7b89897b --- /dev/null +++ b/src/Determinants/usefull_for_ovb.irp.f @@ -0,0 +1,283 @@ + +integer function n_open_shell(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer :: i + n_open_shell = 0 + do i=1,Nint + n_open_shell += popcnt(iand(xor(det_in(i,1),det_in(i,2)),det_in(i,1))) + enddo +end + +integer function n_closed_shell(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer :: i + n_closed_shell = 0 + do i=1,Nint + n_closed_shell += popcnt(iand(det_in(i,1),det_in(i,2))) + enddo +end + +integer function n_closed_shell_cas(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer(bit_kind) :: det_tmp(nint,2) + integer :: i + n_closed_shell_cas = 0 + do i=1,Nint + det_tmp(i,1) = xor(det_in(i,1),reunion_of_core_inact_bitmask(i,1)) + det_tmp(i,2) = xor(det_in(i,2),reunion_of_core_inact_bitmask(i,2)) + enddo +!call debug_det(det_tmp,nint) + do i=1,Nint + n_closed_shell_cas += popcnt(iand(det_tmp(i,1),det_tmp(i,2))) + enddo +end + +subroutine doubly_occ_empty_in_couple(det_in,n_couples,couples,couples_out) + implicit none + use bitmasks + integer, intent(in) :: n_couples,couples(n_couples,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: couples_out(0:n_couples) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + couples_out(0) = .False. + do i = 1, n_couples + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 1)then + couples_out(i) = .False. + cycle + endif + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + couples_out(i) = .False. + else + couples_out(i) = .True. + couples_out(0) = .True. + endif + enddo +end + +subroutine give_index_of_doubly_occ_in_active_space(det_in,doubly_occupied_array) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(N_int,2) + logical, intent(out) :: doubly_occupied_array(n_act_orb) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + do i = 1, n_act_orb + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + i_bite = list_act(i) + call set_bit_to_integer(i_bite,det_tmp_bis,N_int) ! act orb + ! det_tmp is zero except for the orbital "ith" active orbital + integer :: i_count,i_bite + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + doubly_occupied_array(i) = .False. + else + doubly_occupied_array(i) = .True. + endif + enddo +end + +subroutine doubly_occ_empty_in_couple_and_no_hund_elsewhere(det_in,n_couple_no_hund,couple_ion,couple_no_hund,is_ok) + implicit none + use bitmasks + integer, intent(in) :: n_couple_no_hund,couple_ion(2),couple_no_hund(n_couple_no_hund,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: is_ok + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + is_ok = .False. + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couple_ion(1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couple_ion(2),det_tmp_bis,N_int) ! second orb + ! det_tmp is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 1)then + is_ok = .False. + return + endif + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + is_ok = .False. + return + else + do i = 1, n_couple_no_hund + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couple_no_hund (i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couple_no_hund (i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp_bis is zero except for the two orbitals of the couple + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 2)then + is_ok = .False. + return + endif + ! test if orbital there are one alpha and one beta + integer :: i_count_alpha,i_count_beta + i_count_alpha = 0 + i_count_beta = 0 + do j = 1, N_int + i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j))) + i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j))) + enddo + if(i_count_alpha==1.and.i_count_beta==1)then + is_ok = .True. + else + is_ok = .False. + return + endif + enddo + is_ok = .True. + endif +end + + +subroutine neutral_no_hund_in_couple(det_in,n_couples,couples,couples_out) + implicit none + use bitmasks + integer, intent(in) :: n_couples,couples(n_couples,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: couples_out(0:n_couples) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + couples_out(0) = .True. + do i = 1, n_couples + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp_bis is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 2)then + couples_out(i) = .False. + cycle + endif + + ! test if orbital there are one alpha and one beta + integer :: i_count_alpha,i_count_beta + i_count_alpha = 0 + i_count_beta = 0 + do j = 1, N_int + i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j))) + i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j))) + enddo + if(i_count_alpha==1.and.i_count_beta==1)then + couples_out(i) = .True. + else + couples_out(i) = .False. + endif + enddo + do i = 1, n_couples + if(.not.couples_out(i))then + couples_out(0) = .False. + endif + enddo +end + +