mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-24 21:41:46 +02:00
Added the Orbital Entanglement plugin, added the list of active orbitals in bitmasks.irp.f.
This commit is contained in:
parent
3d12ee82da
commit
6d30dabc8b
@ -17,12 +17,12 @@ type keyword =
|
|||||||
| Electrons
|
| Electrons
|
||||||
| Mo_basis
|
| Mo_basis
|
||||||
| Nuclei
|
| Nuclei
|
||||||
| Determinants
|
| Hartree_fock
|
||||||
| Integrals_bielec
|
|
||||||
| Pseudo
|
| Pseudo
|
||||||
|
| Integrals_bielec
|
||||||
| Perturbation
|
| Perturbation
|
||||||
| Properties
|
| Properties
|
||||||
| Hartree_fock
|
| Determinants
|
||||||
;;
|
;;
|
||||||
|
|
||||||
|
|
||||||
@ -32,12 +32,12 @@ let keyword_to_string = function
|
|||||||
| Electrons -> "Electrons"
|
| Electrons -> "Electrons"
|
||||||
| Mo_basis -> "MO basis"
|
| Mo_basis -> "MO basis"
|
||||||
| Nuclei -> "Molecule"
|
| Nuclei -> "Molecule"
|
||||||
| Determinants -> "Determinants"
|
| Hartree_fock -> "Hartree_fock"
|
||||||
| Integrals_bielec -> "Integrals_bielec"
|
|
||||||
| Pseudo -> "Pseudo"
|
| Pseudo -> "Pseudo"
|
||||||
|
| Integrals_bielec -> "Integrals_bielec"
|
||||||
| Perturbation -> "Perturbation"
|
| Perturbation -> "Perturbation"
|
||||||
| Properties -> "Properties"
|
| Properties -> "Properties"
|
||||||
| Hartree_fock -> "Hartree_fock"
|
| Determinants -> "Determinants"
|
||||||
;;
|
;;
|
||||||
|
|
||||||
|
|
||||||
@ -86,18 +86,18 @@ let get s =
|
|||||||
f Ao_basis.(read, to_rst)
|
f Ao_basis.(read, to_rst)
|
||||||
| Determinants_by_hand ->
|
| Determinants_by_hand ->
|
||||||
f Determinants_by_hand.(read, to_rst)
|
f Determinants_by_hand.(read, to_rst)
|
||||||
| Determinants ->
|
| Hartree_fock ->
|
||||||
f Determinants.(read, to_rst)
|
f Hartree_fock.(read, to_rst)
|
||||||
| Integrals_bielec ->
|
|
||||||
f Integrals_bielec.(read, to_rst)
|
|
||||||
| Pseudo ->
|
| Pseudo ->
|
||||||
f Pseudo.(read, to_rst)
|
f Pseudo.(read, to_rst)
|
||||||
|
| Integrals_bielec ->
|
||||||
|
f Integrals_bielec.(read, to_rst)
|
||||||
| Perturbation ->
|
| Perturbation ->
|
||||||
f Perturbation.(read, to_rst)
|
f Perturbation.(read, to_rst)
|
||||||
| Properties ->
|
| Properties ->
|
||||||
f Properties.(read, to_rst)
|
f Properties.(read, to_rst)
|
||||||
| Hartree_fock ->
|
| Determinants ->
|
||||||
f Hartree_fock.(read, to_rst)
|
f Determinants.(read, to_rst)
|
||||||
end
|
end
|
||||||
with
|
with
|
||||||
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
|
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
|
||||||
@ -135,12 +135,12 @@ let set str s =
|
|||||||
in
|
in
|
||||||
let open Input in
|
let open Input in
|
||||||
match s with
|
match s with
|
||||||
| Determinants -> write Determinants.(of_rst, write) s
|
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
|
||||||
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
|
|
||||||
| Pseudo -> write Pseudo.(of_rst, write) s
|
| Pseudo -> write Pseudo.(of_rst, write) s
|
||||||
|
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
|
||||||
| Perturbation -> write Perturbation.(of_rst, write) s
|
| Perturbation -> write Perturbation.(of_rst, write) s
|
||||||
| Properties -> write Properties.(of_rst, write) s
|
| Properties -> write Properties.(of_rst, write) s
|
||||||
| Hartree_fock -> write Hartree_fock.(of_rst, write) s
|
| Determinants -> write Determinants.(of_rst, write) s
|
||||||
| Electrons -> write Electrons.(of_rst, write) s
|
| Electrons -> write Electrons.(of_rst, write) s
|
||||||
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
|
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
|
||||||
| Nuclei -> write Nuclei.(of_rst, write) s
|
| Nuclei -> write Nuclei.(of_rst, write) s
|
||||||
@ -188,12 +188,12 @@ let run check_only ezfio_filename =
|
|||||||
Nuclei ;
|
Nuclei ;
|
||||||
Ao_basis;
|
Ao_basis;
|
||||||
Electrons ;
|
Electrons ;
|
||||||
Determinants ;
|
Hartree_fock ;
|
||||||
Integrals_bielec ;
|
|
||||||
Pseudo ;
|
Pseudo ;
|
||||||
|
Integrals_bielec ;
|
||||||
Perturbation ;
|
Perturbation ;
|
||||||
Properties ;
|
Properties ;
|
||||||
Hartree_fock ;
|
Determinants ;
|
||||||
Mo_basis;
|
Mo_basis;
|
||||||
Determinants_by_hand ;
|
Determinants_by_hand ;
|
||||||
]
|
]
|
||||||
|
1
plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES
Normal file
1
plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES
Normal file
@ -0,0 +1 @@
|
|||||||
|
Determinants
|
345
plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f
Normal file
345
plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f
Normal file
@ -0,0 +1,345 @@
|
|||||||
|
use bitmasks
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [integer, mo_inp_num]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! This is the number of orbitals involved in the entanglement calculation.
|
||||||
|
! It is taken equal to the number of active orbitals n_act_orb.
|
||||||
|
END_DOC
|
||||||
|
mo_inp_num = n_act_orb
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [integer, mo_inp_list, (N_int*bit_kind_size)]
|
||||||
|
&BEGIN_PROVIDER [integer, mo_inp_list_rev, (mo_tot_num)]
|
||||||
|
&BEGIN_PROVIDER [integer(bit_kind), mo_inp_bit_list, (N_int)]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! mo_inp_list is the list of the orbitals involved in the entanglement calculation.
|
||||||
|
! It is taken equal to the list of active orbitals list_act.
|
||||||
|
! mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
|
||||||
|
END_DOC
|
||||||
|
integer :: i
|
||||||
|
|
||||||
|
do i = 1, mo_inp_num
|
||||||
|
mo_inp_list(i)=list_act(i)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do i = 1, mo_inp_num
|
||||||
|
mo_inp_list_rev(mo_inp_list(i))=i
|
||||||
|
enddo
|
||||||
|
call list_to_bitstring( mo_inp_bit_list, mo_inp_list, mo_inp_num, N_int)
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, entropy_one_orb, (mo_inp_num)]
|
||||||
|
&BEGIN_PROVIDER [double precision, entropy_two_orb, (mo_inp_num,mo_inp_num)]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
|
||||||
|
! entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
double precision, allocatable :: ro1(:,:),ro2(:,:,:)
|
||||||
|
integer :: i,j,k,l,ii,jj,iii,istate,kl,info
|
||||||
|
integer, allocatable :: occ(:,:)
|
||||||
|
integer :: n_occ_alpha, n_occ_beta
|
||||||
|
logical, allocatable :: zocc(:,:)
|
||||||
|
logical :: zalpha, zbeta, zalpha2, zbeta2
|
||||||
|
integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,spin1,spin2
|
||||||
|
double precision :: phase
|
||||||
|
integer(bit_kind) :: key_tmp(N_int), key_tmp2(N_int)
|
||||||
|
integer :: ip
|
||||||
|
double precision, parameter :: eps=10.d0**(-14)
|
||||||
|
double precision :: w(16), work(3*16)
|
||||||
|
|
||||||
|
|
||||||
|
allocate(ro1(4,mo_inp_num),ro2(16,16,(mo_inp_num*(mo_inp_num-1)/2)))
|
||||||
|
|
||||||
|
entropy_one_orb = 0.d0
|
||||||
|
entropy_two_orb = 0.d0
|
||||||
|
ro1 = 0.d0
|
||||||
|
ro2 = 0.d0
|
||||||
|
|
||||||
|
allocate (occ(N_int*bit_kind_size,2))
|
||||||
|
allocate (zocc(mo_tot_num,2))
|
||||||
|
|
||||||
|
istate = 1 !Only GS, to be generalized...
|
||||||
|
do ii=1,N_det
|
||||||
|
! We get the occupation of the alpha electrons in occ(:,1)
|
||||||
|
call bitstring_to_list(psi_det(1,1,ii), occ(1,1), n_occ_alpha, N_int)
|
||||||
|
! We get the occupation of the beta electrons in occ(:,2)
|
||||||
|
call bitstring_to_list(psi_det(1,2,ii), occ(1,2), n_occ_beta, N_int)
|
||||||
|
zocc = .false.
|
||||||
|
do i=1,n_occ_alpha
|
||||||
|
zocc(occ(i,1),1)=.true.
|
||||||
|
enddo
|
||||||
|
do i=1,n_occ_beta
|
||||||
|
zocc(occ(i,2),2)=.true.
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do k=1,mo_inp_num
|
||||||
|
zalpha = zocc(mo_inp_list(k),1)
|
||||||
|
zbeta = zocc(mo_inp_list(k),2)
|
||||||
|
! mono start
|
||||||
|
if (zbeta.and.zalpha) then
|
||||||
|
ro1(4,k) = ro1(4,k) + psi_coef(ii,istate)**2 ! double occupied
|
||||||
|
elseif (zalpha) then
|
||||||
|
ro1(2,k) = ro1(2,k) + psi_coef(ii,istate)**2 ! single alpha
|
||||||
|
elseif (zbeta) then
|
||||||
|
ro1(3,k) = ro1(3,k) + psi_coef(ii,istate)**2 ! single beta
|
||||||
|
else
|
||||||
|
ro1(1,k) = ro1(1,k) + psi_coef(ii,istate)**2 ! empty
|
||||||
|
endif
|
||||||
|
! mono stop
|
||||||
|
! double start
|
||||||
|
if (k.eq.mo_inp_num) cycle
|
||||||
|
do l=k+1,mo_inp_num
|
||||||
|
kl=(l-1)*(l-2)/2+k
|
||||||
|
zalpha2 = zocc(mo_inp_list(l),1)
|
||||||
|
zbeta2 = zocc(mo_inp_list(l),2)
|
||||||
|
|
||||||
|
if (zbeta.and.zalpha.and.zbeta2.and.zalpha2) then
|
||||||
|
ro2(16,16,kl) = ro2(16,16,kl) + psi_coef(ii,istate)**2 ! both double occupied
|
||||||
|
else if (zbeta.and.zalpha.and.zbeta2) then
|
||||||
|
ro2(15,15,kl) = ro2(15,15,kl) + psi_coef(ii,istate)**2 ! one double, one beta
|
||||||
|
else if (zbeta.and.zalpha.and.zalpha2) then
|
||||||
|
ro2(13,13,kl) = ro2(13,13,kl) + psi_coef(ii,istate)**2 ! one double, one alpha
|
||||||
|
else if (zbeta.and.zbeta2.and.zalpha2) then
|
||||||
|
ro2(14,14,kl) = ro2(14,14,kl) + psi_coef(ii,istate)**2 ! one beta, one double
|
||||||
|
else if (zalpha.and.zbeta2.and.zalpha2) then
|
||||||
|
ro2(12,12,kl) = ro2(12,12,kl) + psi_coef(ii,istate)**2 ! one alpha, one double
|
||||||
|
else if (zalpha.and.zbeta) then
|
||||||
|
ro2(11,11,kl) = ro2(11,11,kl) + psi_coef(ii,istate)**2 ! one double, one empty
|
||||||
|
else if (zbeta2.and.zalpha2) then
|
||||||
|
ro2(8,8,kl) = ro2(8,8,kl) + psi_coef(ii,istate)**2 ! one empty, one double
|
||||||
|
else if (zbeta.and.zalpha2) then
|
||||||
|
ro2(10,10,kl) = ro2(10,10,kl) + psi_coef(ii,istate)**2 ! one beta, one alpha
|
||||||
|
else if (zalpha.and.zbeta2) then
|
||||||
|
ro2(9,9,kl) = ro2(9,9,kl) + psi_coef(ii,istate)**2 ! one alpha, one beta
|
||||||
|
else if (zbeta.and.zbeta2) then
|
||||||
|
ro2(7,7,kl) = ro2(7,7,kl) + psi_coef(ii,istate)**2 ! one beta, one beta
|
||||||
|
else if (zalpha.and.zalpha2) then
|
||||||
|
ro2(6,6,kl) = ro2(6,6,kl) + psi_coef(ii,istate)**2 ! one alpha, one alpha
|
||||||
|
else if (zbeta) then
|
||||||
|
ro2(5,5,kl) = ro2(5,5,kl) + psi_coef(ii,istate)**2 ! one beta, one empty
|
||||||
|
else if (zbeta2) then
|
||||||
|
ro2(4,4,kl) = ro2(4,4,kl) + psi_coef(ii,istate)**2 ! one empty, one beta
|
||||||
|
else if (zalpha) then
|
||||||
|
ro2(3,3,kl) = ro2(3,3,kl) + psi_coef(ii,istate)**2 ! one alpha, one empty
|
||||||
|
else if (zalpha2) then
|
||||||
|
ro2(2,2,kl) = ro2(2,2,kl) + psi_coef(ii,istate)**2 ! one empty, one alpha
|
||||||
|
else
|
||||||
|
ro2(1,1,kl) = ro2(1,1,kl) + psi_coef(ii,istate)**2 ! both empty
|
||||||
|
end if
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
! stop double
|
||||||
|
|
||||||
|
if (ii.eq.N_det) cycle
|
||||||
|
!Off Diagonal Elements
|
||||||
|
do jj=ii+1,N_det
|
||||||
|
|
||||||
|
call get_excitation_degree(psi_det(1,1,ii),psi_det(1,1,jj),degree,N_int)
|
||||||
|
if (degree.gt.2) cycle
|
||||||
|
ip=0
|
||||||
|
do iii =1,N_int
|
||||||
|
key_tmp(iii) = ior(xor(psi_det(iii,1,ii),psi_det(iii,1,jj)),xor(psi_det(iii,2,ii),psi_det(iii,2,jj)))
|
||||||
|
ip += popcnt(key_tmp(iii))
|
||||||
|
enddo
|
||||||
|
if (ip.ne.2) cycle !They involve more than 2 orbitals.
|
||||||
|
ip=0
|
||||||
|
do iii=1,N_int
|
||||||
|
ip += popcnt(iand(key_tmp(iii),mo_inp_bit_list(iii)))
|
||||||
|
enddo
|
||||||
|
if (ip.ne.2) cycle !They do not involve orbitals of the list.
|
||||||
|
|
||||||
|
if (degree.eq.2) then
|
||||||
|
call get_double_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int)
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2)
|
||||||
|
k=mo_inp_list_rev(h1)
|
||||||
|
l=mo_inp_list_rev(p1)
|
||||||
|
if (k.gt.l) then
|
||||||
|
kl=(k-1)*(k-2)/2+l
|
||||||
|
else
|
||||||
|
kl=(l-1)*(l-2)/2+k
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((.not.zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))&
|
||||||
|
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
|
||||||
|
ro2(8,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(11,8,kl) = ro2(8,11,kl)
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))&
|
||||||
|
.and.(.not.zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
|
||||||
|
ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(10,9,kl) = ro2(9,10,kl)
|
||||||
|
endif
|
||||||
|
if ((zocc(mo_inp_list(k),1)).and.(.not.zocc(mo_inp_list(k),2))&
|
||||||
|
.and.(.not.zocc(mo_inp_list(l),1)).and.(zocc(mo_inp_list(l),2))) then
|
||||||
|
ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(10,9,kl) = ro2(9,10,kl)
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (degree.eq.1) then
|
||||||
|
call get_mono_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int)
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2)
|
||||||
|
k=mo_inp_list_rev(h1)
|
||||||
|
l=mo_inp_list_rev(p1)
|
||||||
|
if (k.gt.l) then
|
||||||
|
kl=(k-1)*(k-2)/2+l
|
||||||
|
else
|
||||||
|
kl=(l-1)*(l-2)/2+k
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((.not.(zocc(mo_inp_list(l),2))).and.&
|
||||||
|
(.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),1))&
|
||||||
|
.and.(.not.(zocc(mo_inp_list(k),2)))) then
|
||||||
|
ro2(2,3,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(3,2,kl)=ro2(2,3,kl)
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((.not.(zocc(mo_inp_list(l),2))).and.&
|
||||||
|
(.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),2))&
|
||||||
|
.and.(.not.(zocc(mo_inp_list(k),1)))) then
|
||||||
|
ro2(4,5,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(5,4,kl)=ro2(4,5,kl)
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
if ((.not.(zocc(mo_inp_list(l),1))).and.& !k doubly occupied, l empty
|
||||||
|
(.not.(zocc(mo_inp_list(l),2))).and.&
|
||||||
|
(zocc(mo_inp_list(k),1)).and.&
|
||||||
|
(zocc(mo_inp_list(k),2))) then
|
||||||
|
if (k.gt.l) then
|
||||||
|
if (spin1.eq.1) then !spin alpha
|
||||||
|
ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(9,8,kl)=ro2(8,9,kl)
|
||||||
|
else !spin beta
|
||||||
|
ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(10,8,kl)=ro2(8,10,kl)
|
||||||
|
endif
|
||||||
|
else ! k.lt.l
|
||||||
|
if (spin1.eq.1) then !spin alpha
|
||||||
|
ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(11,10,kl)=ro2(10,11,kl)
|
||||||
|
else !spin beta
|
||||||
|
ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(11,9,kl)=ro2(9,11,kl)
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((.not.(zocc(mo_inp_list(l),1))).and.& !k alpha, l beta
|
||||||
|
(.not.(zocc(mo_inp_list(k),2))).and.&
|
||||||
|
(zocc(mo_inp_list(k),1)).and.&
|
||||||
|
(zocc(mo_inp_list(l),2))) then
|
||||||
|
if (k.gt.l) then
|
||||||
|
if (spin1.eq.1) then !spin alpha
|
||||||
|
ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(11,10,kl)=ro2(10,11,kl)
|
||||||
|
else !spin beta
|
||||||
|
print*, "problem in k alpha l beta k.gt.l spin beta"
|
||||||
|
endif
|
||||||
|
else ! k.lt.l
|
||||||
|
if (spin1.eq.1) then !spin alpha
|
||||||
|
ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(9,8,kl)=ro2(8,9,kl)
|
||||||
|
else !spin beta
|
||||||
|
print*, "problem in k alpha l beta k.lt.l spin beta"
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ((.not.(zocc(mo_inp_list(k),1))).and.& !k beta, l alpha
|
||||||
|
(.not.(zocc(mo_inp_list(l),2))).and.&
|
||||||
|
(zocc(mo_inp_list(l),1)).and.&
|
||||||
|
(zocc(mo_inp_list(k),2))) then
|
||||||
|
if (k.gt.l) then
|
||||||
|
if (spin1.eq.2) then !spin beta
|
||||||
|
ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(11,9,kl)=ro2(9,11,kl)
|
||||||
|
else !spin alpha
|
||||||
|
print*, "problem in k beta l alpha k.gt.l spin alpha"
|
||||||
|
endif
|
||||||
|
else ! k.lt.l
|
||||||
|
if (spin1.eq.2) then !spin beta
|
||||||
|
ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative
|
||||||
|
ro2(10,8,kl)=ro2(8,10,kl)
|
||||||
|
else !spin alpha
|
||||||
|
print*, "problem in k beta l alpha k.lt.l spin alpha"
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (zocc(mo_inp_list(l),1).and.(.not.zocc(mo_inp_list(l),2))&
|
||||||
|
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
|
||||||
|
ro2(12,13,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(13,12,kl) = ro2(12,13,kl)
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (zocc(mo_inp_list(l),2).and.(.not.zocc(mo_inp_list(l),1))&
|
||||||
|
.and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then
|
||||||
|
ro2(14,15,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate)
|
||||||
|
ro2(15,14,kl) = ro2(14,15,kl)
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
entropy_one_orb=0.d0
|
||||||
|
do k=1,mo_inp_num
|
||||||
|
do i=1,4
|
||||||
|
if (ro1(i,k).ge.eps) then
|
||||||
|
entropy_one_orb(k) = entropy_one_orb(k)-ro1(i,k)*log(ro1(i,k))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
entropy_two_orb=0.d0
|
||||||
|
do k=1,mo_inp_num
|
||||||
|
do l=1,k
|
||||||
|
if (k.eq.l) cycle
|
||||||
|
kl=(k-1)*(k-2)/2+l
|
||||||
|
call dsyev('N','U',16,ro2(1,1,kl),16,w,work,3*16,info)
|
||||||
|
if (info.ne.0) then
|
||||||
|
write(*,*) "Errore in dsyev"
|
||||||
|
endif
|
||||||
|
do j=1,16
|
||||||
|
if (w(j).ge.eps) then
|
||||||
|
entropy_two_orb(k,l) = entropy_two_orb(k,l)-w(j)*log(w(j))
|
||||||
|
entropy_two_orb(l,k) = entropy_two_orb(k,l)
|
||||||
|
elseif ((w(j)).lt.(-eps)) then
|
||||||
|
write(6,*) "Negative Eigenvalue. You have a big problem..."
|
||||||
|
write(6,*) w(j)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
deallocate (occ,zocc,ro1,ro2)
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, mutinf, (mo_inp_num,mo_inp_num)]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
!mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij]
|
||||||
|
!see the refence: 10.1016/j.chemphys.2005.10.018
|
||||||
|
END_DOC
|
||||||
|
integer :: i,j
|
||||||
|
! mutal information:
|
||||||
|
mutinf = 0.d0
|
||||||
|
do i=1,mo_inp_num
|
||||||
|
do j=1,mo_inp_num
|
||||||
|
if (j.eq.i) cycle
|
||||||
|
mutinf(i,j)=-0.5d0*(entropy_two_orb(i,j)-entropy_one_orb(i)-entropy_one_orb(j))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
65
plugins/Orbital_Entanglement/README.rst
Normal file
65
plugins/Orbital_Entanglement/README.rst
Normal file
@ -0,0 +1,65 @@
|
|||||||
|
====================
|
||||||
|
Orbital_Entanglement
|
||||||
|
====================
|
||||||
|
|
||||||
|
Needed Modules
|
||||||
|
==============
|
||||||
|
.. Do not edit this section It was auto-generated
|
||||||
|
.. by the `update_README.py` script.
|
||||||
|
|
||||||
|
|
||||||
|
.. image:: tree_dependency.png
|
||||||
|
|
||||||
|
* `Determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants>`_
|
||||||
|
|
||||||
|
Documentation
|
||||||
|
=============
|
||||||
|
.. Do not edit this section It was auto-generated
|
||||||
|
.. by the `update_README.py` script.
|
||||||
|
|
||||||
|
|
||||||
|
`entropy_one_orb <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L35>`_
|
||||||
|
entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
|
||||||
|
entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
|
||||||
|
|
||||||
|
|
||||||
|
`entropy_two_orb <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L36>`_
|
||||||
|
entropy_one_orb is the one-orbital von Neumann entropy S(1)_i
|
||||||
|
entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij.
|
||||||
|
|
||||||
|
|
||||||
|
`mo_inp_bit_list <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L15>`_
|
||||||
|
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
|
||||||
|
It is taken equal to the list of active orbitals list_act.
|
||||||
|
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
|
||||||
|
|
||||||
|
|
||||||
|
`mo_inp_list <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L13>`_
|
||||||
|
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
|
||||||
|
It is taken equal to the list of active orbitals list_act.
|
||||||
|
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
|
||||||
|
|
||||||
|
|
||||||
|
`mo_inp_list_rev <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L14>`_
|
||||||
|
mo_inp_list is the list of the orbitals involved in the entanglement calculation.
|
||||||
|
It is taken equal to the list of active orbitals list_act.
|
||||||
|
mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i.
|
||||||
|
|
||||||
|
|
||||||
|
`mo_inp_num <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L3>`_
|
||||||
|
This is the number of orbitals involved in the entanglement calculation.
|
||||||
|
It is taken equal to the number of active orbitals n_act_orb.
|
||||||
|
|
||||||
|
|
||||||
|
`mutinf <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f#L330>`_
|
||||||
|
mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij]
|
||||||
|
see the refence: 10.1016/j.chemphys.2005.10.018
|
||||||
|
|
||||||
|
|
||||||
|
`pouet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/print_entanglement.irp.f#L1>`_
|
||||||
|
Undocumented
|
||||||
|
|
||||||
|
|
||||||
|
`routine <http://github.com/LCPQ/quantum_package/tree/master/plugins/Orbital_Entanglement/print_entanglement.irp.f#L9>`_
|
||||||
|
Undocumented
|
||||||
|
|
46
plugins/Orbital_Entanglement/print_entanglement.irp.f
Normal file
46
plugins/Orbital_Entanglement/print_entanglement.irp.f
Normal file
@ -0,0 +1,46 @@
|
|||||||
|
program pouet
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
read_wf = .true.
|
||||||
|
touch read_wf
|
||||||
|
call routine
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine routine
|
||||||
|
implicit none
|
||||||
|
integer:: i,j
|
||||||
|
|
||||||
|
write(6,*) 'Total orbitals: ', mo_tot_num
|
||||||
|
write(6,*) 'Orbitals for entanglement calculation: ', mo_inp_num
|
||||||
|
write(6,*) 'Index: ',(mo_inp_list(i),i=1,mo_inp_num)
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "s1: One-Orbital von Neumann entropy"
|
||||||
|
write(6,'(1000f8.5)') entropy_one_orb
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "s2: Two-Orbital von Neumann entropy"
|
||||||
|
do i=1,mo_inp_num
|
||||||
|
write(6,'(1000f8.5)') (entropy_two_orb(i,j),j=1,mo_inp_num)
|
||||||
|
enddo
|
||||||
|
write(6,*)
|
||||||
|
|
||||||
|
! mutal information:
|
||||||
|
write(6,*) "Mutual Information (Entanglement)"
|
||||||
|
do i=1,mo_inp_num
|
||||||
|
write(6,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
open(17,file=(trim(ezfio_filename)//".entanglement"),status='unknown',form='formatted')
|
||||||
|
|
||||||
|
write(17,'(1000f8.5)') entropy_one_orb
|
||||||
|
do i=1,mo_inp_num
|
||||||
|
write(17,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
close(17)
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "The .entanglement file is the input for the entanglement.py script."
|
||||||
|
write(6,*) "You can find the script in the directory Scripts of QP."
|
||||||
|
write(6,*)
|
||||||
|
end
|
21
plugins/Orbital_Entanglement/tree_dependency
Normal file
21
plugins/Orbital_Entanglement/tree_dependency
Normal file
@ -0,0 +1,21 @@
|
|||||||
|
// ['Orbital_Entanglement']
|
||||||
|
digraph {
|
||||||
|
Orbital_Entanglement [fontcolor=red]
|
||||||
|
Orbital_Entanglement -> Determinants
|
||||||
|
Determinants -> Integrals_Monoelec
|
||||||
|
Integrals_Monoelec -> MO_Basis
|
||||||
|
MO_Basis -> AO_Basis
|
||||||
|
AO_Basis -> Nuclei
|
||||||
|
Nuclei -> Ezfio_files
|
||||||
|
Nuclei -> Utils
|
||||||
|
MO_Basis -> Electrons
|
||||||
|
Electrons -> Ezfio_files
|
||||||
|
Integrals_Monoelec -> Pseudo
|
||||||
|
Pseudo -> Nuclei
|
||||||
|
Determinants -> Integrals_Bielec
|
||||||
|
Integrals_Bielec -> Pseudo
|
||||||
|
Integrals_Bielec -> Bitmask
|
||||||
|
Bitmask -> MO_Basis
|
||||||
|
Integrals_Bielec -> ZMQ
|
||||||
|
ZMQ -> Utils
|
||||||
|
}
|
0
plugins/Orbital_Entanglement/tree_dependency.png
Normal file
0
plugins/Orbital_Entanglement/tree_dependency.png
Normal file
140
scripts/entanglement.py
Executable file
140
scripts/entanglement.py
Executable file
@ -0,0 +1,140 @@
|
|||||||
|
#!/usr/bin/env python
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
import sys
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
from matplotlib import lines
|
||||||
|
from pylab import *
|
||||||
|
|
||||||
|
if __name__ == '__main__':
|
||||||
|
inputfile = sys.argv[1]
|
||||||
|
with open(str(inputfile),'r') as f:
|
||||||
|
vec_s1 = [float(x) for x in f.readline().split()]
|
||||||
|
mat_I = np.array([[float(x) for x in ln.split()] for ln in f])
|
||||||
|
#mat_I = np.loadtxt(str(inputfile),skiprows=1)
|
||||||
|
|
||||||
|
print vec_s1
|
||||||
|
print mat_I
|
||||||
|
#enable the following lines if you dont have a working X display
|
||||||
|
#import matplotlib
|
||||||
|
#matplotlib.use('Agg')
|
||||||
|
|
||||||
|
#MUTUALI INFORMATION PLOT:
|
||||||
|
|
||||||
|
# edit the following line to select the orbitals...
|
||||||
|
zselect=False
|
||||||
|
if zselect:
|
||||||
|
select =[]
|
||||||
|
select.extend(range(34,132))
|
||||||
|
select = [x-1 for x in select]
|
||||||
|
nselect = len(select)
|
||||||
|
vec_s1_s = np.zeros(nselect)
|
||||||
|
mat_I_s = np.zeros((nselect,nselect))
|
||||||
|
print mat_I_s, vec_s1_s
|
||||||
|
|
||||||
|
for i,x in enumerate(select):
|
||||||
|
vec_s1_s[i]=vec_s1[x]
|
||||||
|
for i,x in enumerate(select):
|
||||||
|
for j,y in enumerate(select):
|
||||||
|
mat_I_s[i,j] = mat_I[x,y]
|
||||||
|
|
||||||
|
|
||||||
|
print vec_s1_s
|
||||||
|
print mat_I_s
|
||||||
|
vec_s1 = vec_s1_s
|
||||||
|
mat_I = mat_I_s
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#end selection
|
||||||
|
|
||||||
|
|
||||||
|
N= len(mat_I)
|
||||||
|
print N
|
||||||
|
theta = np.zeros(N)
|
||||||
|
r=np.zeros(N)
|
||||||
|
labels=np.zeros(N)
|
||||||
|
area=np.zeros(N)
|
||||||
|
|
||||||
|
for i in range(N):
|
||||||
|
theta[i]=-2*pi/N*i+pi/2
|
||||||
|
r[i]=1.0
|
||||||
|
if zselect:
|
||||||
|
labels[i]=select[i]+1
|
||||||
|
else:
|
||||||
|
labels[i]=i+1
|
||||||
|
area[i]=vec_s1[i]*500
|
||||||
|
|
||||||
|
if (len(sys.argv) > 2):
|
||||||
|
for j in range(2,len(sys.argv)):
|
||||||
|
if (sys.argv[j] == '-c'):
|
||||||
|
coordfile = sys.argv[j+1]
|
||||||
|
with open(str(coordfile),'r') as f:
|
||||||
|
for i in range(N):
|
||||||
|
line = f.readline().split()
|
||||||
|
theta[i] = float(line[0]) +pi/2
|
||||||
|
r[i] = float(line[1])
|
||||||
|
print 'theta=',theta
|
||||||
|
print 'r=',r
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
ax = plt.subplot(111,polar=True)
|
||||||
|
ax.set_xticklabels([])
|
||||||
|
ax.set_yticklabels([])
|
||||||
|
ax.grid(b=False)
|
||||||
|
c = plt.scatter(theta,r,c="Red",s=area)
|
||||||
|
|
||||||
|
|
||||||
|
if (len(sys.argv) > 2):
|
||||||
|
for j in range(2,len(sys.argv)):
|
||||||
|
print 'sys.argv',j,sys.argv[j]
|
||||||
|
if (sys.argv[j] == '-t'):
|
||||||
|
plt.title(sys.argv[j+1])
|
||||||
|
break
|
||||||
|
else:
|
||||||
|
plt.title('N = '+str(N/2)+', Results file: '+sys.argv[1])
|
||||||
|
|
||||||
|
#this is dummy:
|
||||||
|
#c1 = plt.scatter(theta,(r+0.05),c="red",s=0)
|
||||||
|
|
||||||
|
for i in range(N):
|
||||||
|
# plt.annotate(int(labels[i]),xy=(theta[i],(r[i]+0.2)),size='xx-large',)
|
||||||
|
plt.text(theta[i],(r[i]+0.25),int(labels[i]),size='x-small',ha='center',va='center')
|
||||||
|
for j in range(i,N):
|
||||||
|
x =[theta[i],theta[j]]
|
||||||
|
y =[r[i],r[j]]
|
||||||
|
#y =[1,1]
|
||||||
|
if mat_I[i,j] >= 0.1:
|
||||||
|
line1 = lines.Line2D(x, y, linewidth=3, color='black',linestyle='-', alpha=0.8,label='0.1')
|
||||||
|
ax.add_line(line1)
|
||||||
|
elif mat_I[i,j] >=0.01:
|
||||||
|
line2 = lines.Line2D(x, y, linewidth=3, color='green',linestyle='-', alpha=0.8,label='0.01')
|
||||||
|
ax.add_line(line2)
|
||||||
|
elif mat_I[i,j] >=0.001:
|
||||||
|
line3 = lines.Line2D(x, y, linewidth=3, color='gray',linestyle='-', alpha=0.6,label='0.001')
|
||||||
|
ax.add_line(line3)
|
||||||
|
params = {'legend.fontsize': 16,
|
||||||
|
'legend.linewidth': 2}
|
||||||
|
|
||||||
|
#need to check if the three types of lines really exist.
|
||||||
|
if 'line1' in locals() and 'line2' in locals() and 'line3' in locals():
|
||||||
|
ax.legend([line1, line2, line3],[line1.get_label(),line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line1' in locals() and 'line2' in locals():
|
||||||
|
ax.legend([line1, line2],[line1.get_label(),line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line1' in locals() and 'line3' in locals():
|
||||||
|
ax.legend([line1, line3],[line1.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line2' in locals() and 'line3' in locals():
|
||||||
|
ax.legend([line2, line3],[line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line1' in locals():
|
||||||
|
ax.legend([line1],[line1.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line2' in locals():
|
||||||
|
ax.legend([line3],[line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
elif 'line3' in locals():
|
||||||
|
ax.legend([line2],[line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True)
|
||||||
|
#probably not the best way to code it.
|
||||||
|
|
||||||
|
#for saving, enable matplotlib.use('Agg')
|
||||||
|
#plt.savefig(str(sys.argv[1])+'.eps')
|
||||||
|
plt.show()
|
@ -407,3 +407,31 @@ END_PROVIDER
|
|||||||
unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2))
|
unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2))
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ integer, n_act_orb]
|
||||||
|
BEGIN_DOC
|
||||||
|
! number of active orbitals
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i,j
|
||||||
|
n_act_orb = 0
|
||||||
|
do i = 1, N_int
|
||||||
|
n_act_orb += popcnt(cas_bitmask(i,1,1))
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [integer, list_act, (n_act_orb)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! list of active orbitals
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: occ_act(N_int*bit_kind_size)
|
||||||
|
integer :: itest,i
|
||||||
|
occ_act = 0
|
||||||
|
call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int)
|
||||||
|
ASSERT(itest==n_act_orb)
|
||||||
|
do i = 1, n_act_orb
|
||||||
|
list_act(i) = occ_act(i)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
Loading…
Reference in New Issue
Block a user