mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-22 20:35:19 +01:00
511 lines
19 KiB
Fortran
511 lines
19 KiB
Fortran
|
|
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
|
|
|