LCPQ
/
quantum_package
mirror of https://github.com/LCPQ/quantum_package
10
0
Fork 0
Code Issues Releases Wiki Activity

added print_hcc and OVB plugin

This commit is contained in:
Emmanuel Giner 2016-02-17 10:52:57 +01:00
parent 19e276fc0d
commit 27a121bc1b
12 changed files with 967 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -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

7
plugins/Hartree_Fock/EZFIO.cfg
View File

@ -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

2
plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Integrals_Bielec MOGuess
Integrals_Bielec MOGuess Bitmask

74
plugins/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -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

1
plugins/OVB/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Determinants Psiref_CAS

20
plugins/OVB/README.rst Normal file
View File

@ -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.

510
plugins/OVB/ovb_components.irp.f Normal file
View File

@ -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

27
plugins/OVB/print_ovb.irp.f Normal file
View File

@ -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

17
plugins/Properties/print_hcc.irp.f Normal file
View File

@ -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

21
plugins/Psiref_Utils/psi_ref_utils.irp.f
View File

@ -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

61
src/Bitmask/bitmasks.irp.f
View File

@ -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 ]

283
src/Determinants/usefull_for_ovb.irp.f Normal file
View File

@ -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