From 6d30dabc8b6aa7cbd983a6897f743f7e1e5352b2 Mon Sep 17 00:00:00 2001 From: Lorenzo Tenti Date: Fri, 12 Feb 2016 12:27:13 +0100 Subject: [PATCH] Added the Orbital Entanglement plugin, added the list of active orbitals in bitmasks.irp.f. --- ocaml/qp_edit.ml | 36 +- .../NEEDED_CHILDREN_MODULES | 1 + .../Orbital_Entanglement.irp.f | 345 ++++++++++++++++++ plugins/Orbital_Entanglement/README.rst | 65 ++++ .../print_entanglement.irp.f | 46 +++ plugins/Orbital_Entanglement/tree_dependency | 21 ++ .../Orbital_Entanglement/tree_dependency.png | 0 scripts/entanglement.py | 140 +++++++ src/Bitmask/bitmasks.irp.f | 28 ++ 9 files changed, 664 insertions(+), 18 deletions(-) create mode 100644 plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f create mode 100644 plugins/Orbital_Entanglement/README.rst create mode 100644 plugins/Orbital_Entanglement/print_entanglement.irp.f create mode 100644 plugins/Orbital_Entanglement/tree_dependency create mode 100644 plugins/Orbital_Entanglement/tree_dependency.png create mode 100755 scripts/entanglement.py diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml index a693aa2f..53e3ea59 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -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 ; ] diff --git a/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES b/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..aae89501 --- /dev/null +++ b/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants diff --git a/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f b/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f new file mode 100644 index 00000000..74d7d8ae --- /dev/null +++ b/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f @@ -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 diff --git a/plugins/Orbital_Entanglement/README.rst b/plugins/Orbital_Entanglement/README.rst new file mode 100644 index 00000000..9b5b8119 --- /dev/null +++ b/plugins/Orbital_Entanglement/README.rst @@ -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 `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`entropy_one_orb `_ + 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 `_ + 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 `_ + 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 `_ + 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 `_ + 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 `_ + 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 `_ + 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 `_ + Undocumented + + +`routine `_ + Undocumented + diff --git a/plugins/Orbital_Entanglement/print_entanglement.irp.f b/plugins/Orbital_Entanglement/print_entanglement.irp.f new file mode 100644 index 00000000..e7e732cc --- /dev/null +++ b/plugins/Orbital_Entanglement/print_entanglement.irp.f @@ -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 diff --git a/plugins/Orbital_Entanglement/tree_dependency b/plugins/Orbital_Entanglement/tree_dependency new file mode 100644 index 00000000..a601745c --- /dev/null +++ b/plugins/Orbital_Entanglement/tree_dependency @@ -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 +} \ No newline at end of file diff --git a/plugins/Orbital_Entanglement/tree_dependency.png b/plugins/Orbital_Entanglement/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/scripts/entanglement.py b/scripts/entanglement.py new file mode 100755 index 00000000..c55cc98e --- /dev/null +++ b/scripts/entanglement.py @@ -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() diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index f54532ca..70a84a42 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -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