mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Merge pull request #141 from lorenzotenti/master
Added the Orbital Entanglement plugin
This commit is contained in:
Thomas Applencourt
commit
023ff8cefb
@ -17,12 +17,12 @@ type keyword =
|
||||
| Electrons
|
||||
| Mo_basis
|
||||
| Nuclei
|
||||
| Determinants
|
||||
| Integrals_bielec
|
||||
| Hartree_fock
|
||||
| Pseudo
|
||||
| Integrals_bielec
|
||||
| Perturbation
|
||||
| Properties
|
||||
| Hartree_fock
|
||||
| Determinants
|
||||
;;
|
||||
|
||||
|
||||
@ -32,12 +32,12 @@ let keyword_to_string = function
|
||||
| Electrons -> "Electrons"
|
||||
| Mo_basis -> "MO basis"
|
||||
| Nuclei -> "Molecule"
|
||||
| Determinants -> "Determinants"
|
||||
| Integrals_bielec -> "Integrals_bielec"
|
||||
| Hartree_fock -> "Hartree_fock"
|
||||
| Pseudo -> "Pseudo"
|
||||
| Integrals_bielec -> "Integrals_bielec"
|
||||
| Perturbation -> "Perturbation"
|
||||
| Properties -> "Properties"
|
||||
| Hartree_fock -> "Hartree_fock"
|
||||
| Determinants -> "Determinants"
|
||||
;;
|
||||
|
||||
|
||||
@ -86,18 +86,18 @@ let get s =
|
||||
f Ao_basis.(read, to_rst)
|
||||
| Determinants_by_hand ->
|
||||
f Determinants_by_hand.(read, to_rst)
|
||||
| Determinants ->
|
||||
f Determinants.(read, to_rst)
|
||||
| Integrals_bielec ->
|
||||
f Integrals_bielec.(read, to_rst)
|
||||
| Hartree_fock ->
|
||||
f Hartree_fock.(read, to_rst)
|
||||
| Pseudo ->
|
||||
f Pseudo.(read, to_rst)
|
||||
| Integrals_bielec ->
|
||||
f Integrals_bielec.(read, to_rst)
|
||||
| Perturbation ->
|
||||
f Perturbation.(read, to_rst)
|
||||
| Properties ->
|
||||
f Properties.(read, to_rst)
|
||||
| Hartree_fock ->
|
||||
f Hartree_fock.(read, to_rst)
|
||||
| Determinants ->
|
||||
f Determinants.(read, to_rst)
|
||||
end
|
||||
with
|
||||
| Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
|
||||
@ -135,12 +135,12 @@ let set str s =
|
||||
in
|
||||
let open Input in
|
||||
match s with
|
||||
| Determinants -> write Determinants.(of_rst, write) s
|
||||
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
|
||||
| Hartree_fock -> write Hartree_fock.(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
|
||||
| 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
|
||||
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
|
||||
| Nuclei -> write Nuclei.(of_rst, write) s
|
||||
@ -188,12 +188,12 @@ let run check_only ezfio_filename =
|
||||
Nuclei ;
|
||||
Ao_basis;
|
||||
Electrons ;
|
||||
Determinants ;
|
||||
Integrals_bielec ;
|
||||
Hartree_fock ;
|
||||
Pseudo ;
|
||||
Integrals_bielec ;
|
||||
Perturbation ;
|
||||
Properties ;
|
||||
Hartree_fock ;
|
||||
Determinants ;
|
||||
Mo_basis;
|
||||
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))
|
||||
enddo
|
||||
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