LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
Code Issues Releases Wiki Activity

Merge pull request #146 from eginer/master

Browse Source 
add the OVB analysis and the mulliken and hyperfine coupling constants
This commit is contained in:
Thomas Applencourt 2016-02-17 12:12:33 +01:00
commit 6b8ef16d0a
23 changed files with 1374 additions and 70 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

2
plugins/Hartree_Fock/.gitignore vendored
View File

@ -5,7 +5,6 @@ AO_Basis
Bitmask
Electrons
Ezfio_files
Huckel_guess
IRPF90_man
IRPF90_temp
Integrals_Bielec
@ -16,7 +15,6 @@ Makefile
Makefile.depend
Nuclei
Pseudo
SCF
Utils
ZMQ
ezfio_interface.irp.f

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

135
plugins/Properties/hyperfine_constants.irp.f Normal file
View File

@ -0,0 +1,135 @@
BEGIN_PROVIDER [double precision, spin_density_at_nucleous, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
do k = 1, ao_num
accu += one_body_spin_density_ao(k,j) * aos_array(k) * aos_array(j)
enddo
enddo
spin_density_at_nucleous(i) = accu
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_from_mo, (nucl_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_per_mo, (nucl_num,mo_tot_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: contrib
double precision :: mo_values(mo_tot_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
spin_density_at_nucleous_from_mo(i) = 0.d0
do k = 1, mo_tot_num
mo_values(k) = 0.d0
do j = 1, ao_num
mo_values(k) += mo_coef(j,k) * aos_array(j)
enddo
enddo
do k = 1, mo_tot_num
spin_density_at_nucleous_contrib_per_mo(i,k) = 0.d0
do m = 1, mo_tot_num
contrib = one_body_spin_density_mo(k,m) * mo_values(k) * mo_values(m)
spin_density_at_nucleous_from_mo(i) += contrib
spin_density_at_nucleous_contrib_per_mo(i,k) += contrib
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo, (nucl_num,mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo_test, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
spin_density_at_nucleous_contrib_mo_test = 0.d0
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: c_i1,c_j1
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do k = 1, mo_tot_num
do m = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
c_i1 = mo_coef(j,k)
do l = 1, ao_num
c_j1 = c_i1*mo_coef(l,m)
accu += one_body_spin_density_mo(k,m) * aos_array(l) * aos_array(j) * c_j1
enddo
enddo
spin_density_at_nucleous_contrib_mo(i,k,m) = accu
spin_density_at_nucleous_contrib_mo_test(i) += accu
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, conversion_factor_mhz_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_gauss_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_cm_1_hcc, (100)]
BEGIN_DOC
! Conversion factor for the calculation of the hcc, according to the nuclear charge
END_DOC
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
! hydrogen
conversion_factor_mhz_hcc(1) = 4469.84692227102460d0
conversion_factor_gauss_hcc(1) = 1594.95296390862904d0
conversion_factor_cm_1_hcc(1) = 1490.98044430157870d0
! Li
conversion_factor_mhz_hcc(3) = 1737.2746512855997d0
conversion_factor_gauss_hcc(3) = 619.9027742370165d0
conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
! carbon
conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
conversion_factor_gauss_hcc(6) = 401.136570647523058d0
conversion_factor_cm_1_hcc(6) = 374.987097339830086d0
! nitrogen
conversion_factor_mhz_hcc(7) = 323.102093833793390d0
conversion_factor_gauss_hcc(7) = 115.290892768082614d0
conversion_factor_cm_1_hcc(7) = 107.775257586297698d0
! Oxygen
conversion_factor_mhz_hcc(8) = -606.1958551736545d0
conversion_factor_gauss_hcc(8) = -216.30574771560407d0
conversion_factor_cm_1_hcc(8) = -202.20517197179822d0
END_PROVIDER
BEGIN_PROVIDER [double precision, iso_hcc_mhz, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_gauss, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_cm_1, (nucl_num)]
BEGIN_DOC
! isotropic hyperfine coupling constants among the various atoms
END_DOC
integer :: i
do i = 1, nucl_num
iso_hcc_mhz(i) = conversion_factor_mhz_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !* 0.5d0
iso_hcc_gauss(i) = conversion_factor_gauss_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)!* 0.5d0
iso_hcc_cm_1(i) = conversion_factor_cm_1_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !*0.5d0
enddo
END_PROVIDER

107
plugins/Properties/mulliken.irp.f Normal file
View File

@ -0,0 +1,107 @@
BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
spin_population = 0.d0
do i = 1, ao_num
do j = 1, ao_num
spin_population(j,i) = one_body_spin_density_ao(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)]
implicit none
integer :: i
double precision :: accu
spin_population_angular_momentum = 0.d0
do i = 1, ao_num
spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_gross_orbital_product, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product for the spin population
END_DOC
do i = 1, ao_num
do j = 1, ao_num
spin_gross_orbital_product(i) += spin_population(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_spin_densities, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!ATOMIC SPIN POPULATION (ALPHA MINUS BETA)
END_DOC
mulliken_spin_densities = 0.d0
do i = 1, ao_num
mulliken_spin_densities(ao_nucl(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, electronic_population_alpha, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [double precision, electronic_population_beta, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
electronic_population_alpha = 0.d0
electronic_population_beta = 0.d0
do i = 1, ao_num
do j = 1, ao_num
electronic_population_alpha(j,i) = one_body_dm_ao_alpha(i,j) * ao_overlap(i,j)
electronic_population_beta(j,i) = one_body_dm_ao_beta(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, gross_orbital_product_alpha, (ao_num)]
&BEGIN_PROVIDER [double precision, gross_orbital_product_beta, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product
END_DOC
do i = 1, ao_num
do j = 1, ao_num
gross_orbital_product_alpha(i) += electronic_population_alpha(j,i)
gross_orbital_product_beta(i) += electronic_population_beta(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_densities_alpha, (nucl_num)]
&BEGIN_PROVIDER [double precision, mulliken_densities_beta, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!
END_DOC
mulliken_densities_alpha = 0.d0
mulliken_densities_beta = 0.d0
do i = 1, ao_num
mulliken_densities_alpha(ao_nucl(i)) += gross_orbital_product_alpha(i)
mulliken_densities_beta(ao_nucl(i)) += gross_orbital_product_beta(i)
enddo
END_PROVIDER

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

35
plugins/Properties/print_mulliken.irp.f Normal file
View File

@ -0,0 +1,35 @@
program print_mulliken
implicit none
read_wf = .True.
touch read_wf
print*,'Mulliken spin densities'
call test
end
subroutine test
double precision :: accu
integer :: i
integer :: j
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += spin_gross_orbital_product(i)
write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += spin_population_angular_momentum(i)
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu
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

4
scripts/generate_h_apply.py
View File

@ -156,11 +156,11 @@ class H_apply(object):
def filter_only_1h1p(self):
self["filter_only_1h1p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(hole).eq..False.) cycle
if (is_a_1h1p(hole).eqv..False.) cycle
"""
self["filter_only_1h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(key).eq..False.) cycle
if (is_a_1h1p(key).eqv..False.) cycle
"""

27
src/AO_Basis/aos.irp.f
View File

@ -146,3 +146,30 @@ integer function ao_power_index(nx,ny,nz)
ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
end
BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
&BEGIN_PROVIDER [ integer, ao_l_max ]
&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ]
implicit none
BEGIN_DOC
! ao_l = l value of the AO: a+b+c in x^a y^b z^c
END_DOC
integer :: i
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
ao_l_char(i) = l_to_charater(ao_l(i))
enddo
ao_l_max = maxval(ao_l)
END_PROVIDER
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
BEGIN_DOC
! character corresponding to the "L" value of an AO orbital
END_DOC
implicit none
l_to_charater(0)='S'
l_to_charater(1)='P'
l_to_charater(2)='D'
l_to_charater(3)='F'
l_to_charater(4)='G'
END_PROVIDER

48
src/AO_Basis/aos_value.irp.f Normal file
View File

@ -0,0 +1,48 @@
double precision function ao_value(i,r)
implicit none
BEGIN_DOC
! return the value of the ith ao at point r
END_DOC
double precision, intent(in) :: r(3)
integer, intent(in) :: i
integer :: m,num_ao
double precision :: center_ao(3)
double precision :: beta
integer :: power_ao(3)
num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3)
double precision :: accu,dx,dy,dz,r2
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1)
dy = dy**power_ao(2)
dz = dz**power_ao(3)
accu = 0.d0
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
enddo
ao_value = accu * dx * dy * dz
end
subroutine give_all_aos_at_r(r,aos_array)
implicit none
BEGIN_dOC
! gives the values of aos at a given point r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
integer :: i
double precision :: ao_value
do i = 1, ao_num
aos_array(i) = ao_value(i,r)
enddo
end

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 ]

51
src/Determinants/density_matrix.irp.f
View File

@ -206,3 +206,54 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
state_average_weight = 1.d0/dble(N_states)
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_spin_density_ao, (ao_num_align,ao_num) ]
BEGIN_DOC
! one body spin density matrix on the AO basis : rho_AO(alpha) - rho_AO(beta)
END_DOC
implicit none
integer :: i,j,k,l
double precision :: dm_mo
one_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_spin_density_mo(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_spin_density_ao(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ]
BEGIN_DOC
! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta)
END_DOC
implicit none
integer :: i,j,k,l
double precision :: dm_mo
one_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_dm_mo_alpha(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
enddo
enddo
enddo
enddo
END_PROVIDER

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

3
src/Integrals_Bielec/.gitignore vendored
View File

@ -17,5 +17,4 @@ ZMQ
ezfio_interface.irp.f
irpf90.make
irpf90_entities
tags
test_integrals
tags

2
src/Integrals_Monoelec/.gitignore vendored
View File

@ -12,9 +12,7 @@ Makefile.depend
Nuclei
Pseudo
Utils
check_orthonormality
ezfio_interface.irp.f
irpf90.make
irpf90_entities
save_ortho_mos
tags

5
src/MOGuess/.gitignore vendored
View File

@ -4,7 +4,6 @@
AO_Basis
Electrons
Ezfio_files
H_CORE_guess
IRPF90_man
IRPF90_temp
Integrals_Monoelec
@ -15,8 +14,6 @@ Nuclei
Pseudo
Utils
ezfio_interface.irp.f
guess_overlap
irpf90.make
irpf90_entities
tags
truncate_mos
tags