From 1c42e50cb86d737e0bbf842b8a27107d69a0d42d Mon Sep 17 00:00:00 2001 From: SamueleGiuli Date: Wed, 22 Mar 2023 10:22:04 +0100 Subject: [PATCH 1/4] Small Bugs Fixed in excitonic forcing field --- src/ED_NONSU2/stored/Himp.f90 | 4 ++-- src/ED_NORMAL/ED_HAMILTONIAN_NORMAL_STORED_HxV.f90 | 1 + src/ED_NORMAL/stored/H_dw.f90 | 6 +++--- src/ED_NORMAL/stored/H_up.f90 | 4 ++-- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/ED_NONSU2/stored/Himp.f90 b/src/ED_NONSU2/stored/Himp.f90 index e72a46d..c6cfb1d 100644 --- a/src/ED_NONSU2/stored/Himp.f90 +++ b/src/ED_NONSU2/stored/Himp.f90 @@ -104,7 +104,7 @@ !Evaluate: F.dot.T = sum_a={0,x,y,z} F_a.T^a if(any(exc_field/=0d0))then do iorb=1,Norb - do jorb=iorb+1,Norb + do jorb=1,Norb ! !F_0. T^0_ab := F_0 . (C^+_{a,up}C_{b,up} + C^+_{a,dw}C_{b,dw}) !F_z. T^z_ab := F_z . (C^+_{a,up}C_{b,up} - C^+_{a,dw}C_{b,dw}) @@ -140,7 +140,7 @@ call cdg(iorb,k1,k2,sg2) j = binary_search(Hsector%H(1)%map,k2) ! - htmp = one*exc_field(1)*sg1*sg2 + htmp = exc_field(1)*sg1*sg2 select case(MpiStatus) case (.true.) call sp_insert_element(MpiComm,spH0,htmp,i,j) diff --git a/src/ED_NORMAL/ED_HAMILTONIAN_NORMAL_STORED_HxV.f90 b/src/ED_NORMAL/ED_HAMILTONIAN_NORMAL_STORED_HxV.f90 index 01f5273..aabac80 100644 --- a/src/ED_NORMAL/ED_HAMILTONIAN_NORMAL_STORED_HxV.f90 +++ b/src/ED_NORMAL/ED_HAMILTONIAN_NORMAL_STORED_HxV.f90 @@ -560,6 +560,7 @@ subroutine spMatVec_normal_orbs(Nloc,v,Hv) ! Hv=0d0 ! + do i = 1,Nloc i_el = mod(i-1,DimUp*DimDw) + 1 ! diff --git a/src/ED_NORMAL/stored/H_dw.f90 b/src/ED_NORMAL/stored/H_dw.f90 index e5678cd..7373070 100644 --- a/src/ED_NORMAL/stored/H_dw.f90 +++ b/src/ED_NORMAL/stored/H_dw.f90 @@ -85,14 +85,14 @@ !F_z. T^z_ab := F_z . (- C^+_{a,dw}C_{b,dw}) if(any(exc_field/=0d0))then do iorb=1,Norb - do jorb=iorb+1,Norb + do jorb=1,Norb Jcondition = (Ndw(jorb)==1) .AND. (Ndw(iorb)==0) if (Jcondition) then call c(jorb,mdw,k1,sg1) call cdg(iorb,k1,k2,sg2) - iup = binary_search(Hsector%H(2)%map,k2) + idw = binary_search(Hsector%H(2)%map,k2) ! - htmp = exc_field(1)*sg1*sg2 + htmp = exc_field(1)*sg1*sg2 call sp_insert_element(spH0dws(1),htmp,idw,jdw) ! htmp = -exc_field(4)*sg1*sg2 diff --git a/src/ED_NORMAL/stored/H_up.f90 b/src/ED_NORMAL/stored/H_up.f90 index 22e08c6..980b1ec 100644 --- a/src/ED_NORMAL/stored/H_up.f90 +++ b/src/ED_NORMAL/stored/H_up.f90 @@ -85,14 +85,14 @@ !F_z. T^z_ab := F_z . C^+_{a,up}C_{b,up} if(any(exc_field/=0d0))then do iorb=1,Norb - do jorb=iorb+1,Norb + do jorb=1,Norb Jcondition = (Nup(jorb)==1) .AND. (Nup(iorb)==0) if (Jcondition) then call c(jorb,mup,k1,sg1) call cdg(iorb,k1,k2,sg2) iup = binary_search(Hsector%H(1)%map,k2) ! - htmp = exc_field(1)*sg1*sg2 + htmp = exc_field(1)*sg1*sg2 call sp_insert_element(spH0ups(1),htmp,iup,jup) ! htmp = exc_field(4)*sg1*sg2 From b7f4b8c295d2d349557bf5ac26321db7ae13bfb5 Mon Sep 17 00:00:00 2001 From: SamueleGiuli Date: Tue, 16 May 2023 11:58:50 +0200 Subject: [PATCH 2/4] Minor Update for ED_NORMAL and ED_SUPERC: Added external field "A_ph" coupled to phononic displacement operator (b+bdg)/sqrt(2) Added calculation of and that are now printed inside imp_last and imp_all --- src/ED_INPUT_VARS.f90 | 2 ++ src/ED_NORMAL/ED_OBSERVABLES_NORMAL.f90 | 20 +++++++++++++++++--- src/ED_NORMAL/stored/H_ph.f90 | 14 +++++++++++++- src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 | 19 ++++++++++++++++--- src/ED_SUPERC/stored/H_ph.f90 | 13 +++++++++++++ 5 files changed, 61 insertions(+), 7 deletions(-) diff --git a/src/ED_INPUT_VARS.f90 b/src/ED_INPUT_VARS.f90 index e43cc69..1d65576 100644 --- a/src/ED_INPUT_VARS.f90 +++ b/src/ED_INPUT_VARS.f90 @@ -25,6 +25,7 @@ MODULE ED_INPUT_VARS complex(8),allocatable :: g_ph(:,:) !g_ph: electron-phonon coupling constant all real(8),allocatable :: g_ph_diag(:) !g_ph: electron-phonon coupling constant diagonal (density) real(8) :: w0_ph !w0_ph: phonon frequency (constant) + real(8) :: A_ph !A_ph: phonon field coupled to displacement operator (constant) ! integer :: Nsuccess !# of repeated success to fall below convergence threshold real(8) :: dmft_error !dmft convergence threshold @@ -171,6 +172,7 @@ subroutine ed_read_input(INPUTunit) allocate(g_ph_diag(Norb)) ! THIS SHOULD BE A MATRIX Norb*Norb call parse_input_variable(g_ph_diag,"G_PH",INPUTunit,default=(/( 0d0,i=1,Norb )/),comment="Electron-phonon coupling density constant") call parse_input_variable(w0_ph,"W0_PH",INPUTunit,default=0.d0,comment="Phonon frequency") + call parse_input_variable(A_ph,"A_PH",INPUTunit,default=0.d0,comment="Forcing field coupled to phonon's displacement operator") call parse_input_variable(GPHfile,"GPHfile",INPUTunit,default="NONE",comment="File of Phonon couplings. Put NONE to use only density couplings.") allocate(spin_field_x(Norb)) diff --git a/src/ED_NORMAL/ED_OBSERVABLES_NORMAL.f90 b/src/ED_NORMAL/ED_OBSERVABLES_NORMAL.f90 index 147e1ef..a8c9305 100644 --- a/src/ED_NORMAL/ED_OBSERVABLES_NORMAL.f90 +++ b/src/ED_NORMAL/ED_OBSERVABLES_NORMAL.f90 @@ -26,6 +26,7 @@ MODULE ED_OBSERVABLES_NORMAL real(8),dimension(:,:),allocatable :: exct_tz real(8),dimension(:,:),allocatable :: zimp,simp real(8) :: dens_ph + real(8) :: X_ph, X2_ph real(8) :: s2tot real(8) :: Egs real(8) :: Ei @@ -102,6 +103,8 @@ subroutine observables_normal() Prob = 0.d0 prob_ph = 0.d0 dens_ph = 0.d0 + X_ph = 0.d0 + X2_ph= 0.d0 pdf_ph = 0.d0 pdf_part= 0.d0 w_ph = w0_ph @@ -171,6 +174,17 @@ subroutine observables_normal() i_el = mod(i-1,sectorI%DimEl) + 1 prob_ph(iph) = prob_ph(iph) + gs_weight dens_ph = dens_ph + (iph-1)*gs_weight + + ! and with X=(b+bdg)/sqrt(2) + if(iph1 .AND. iph==1) then @@ -736,7 +750,7 @@ subroutine write_legend() write(unit,"(A1,90(A10,6X))") "#",((reg(txtfy(iorb+(ispin-1)*Norb))//"sig_"//reg(txtfy(iorb))//"s"//reg(txtfy(ispin)),iorb=1,Norb),ispin=1,Nspin) write(unit,"(A1,90(A10,6X))") "# *****" write(unit,"(A1,90(A10,6X))") "# imp_last.ed" - write(unit,"(A1,90(A10,6X))") "#", "1s2tot", "2egs", "3nph", "4w_ph" + write(unit,"(A1,90(A10,6X))") "#", "1s2tot", "2egs", "3nph", "4w_ph","5X_ph", "6X2_ph" write(unit,"(A1,90(A10,6X))") "# *****" write(unit,"(A1,90(A10,6X))") "# exciton_last.ed" write(unit,"(A1,90(A10,6X))") "#","1S_0" , "2T_z" @@ -828,7 +842,7 @@ subroutine write_observables() ! unit = free_unit() open(unit,file="imp_all"//reg(ed_file_suffix)//".ed",position='append') - write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph + write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph, X_ph, X2_ph close(unit) ! unit = free_unit() @@ -888,7 +902,7 @@ subroutine write_observables() ! unit = free_unit() open(unit,file="imp_last"//reg(ed_file_suffix)//".ed") - write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph + write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph, X_ph, X2_ph close(unit) ! unit = free_unit() diff --git a/src/ED_NORMAL/stored/H_ph.f90 b/src/ED_NORMAL/stored/H_ph.f90 index 9fdfff6..844218d 100644 --- a/src/ED_NORMAL/stored/H_ph.f90 +++ b/src/ED_NORMAL/stored/H_ph.f90 @@ -1,6 +1,18 @@ - !Phononic hamiltonian H_ph = w0 b^+ b + !Phononic hamiltonian H_ph = w0 b^+ b + A(bdg + b)/sqrt(2) !Phonon coupled to densities : phi_i * n_{i,iorb,sigma} do iph=1,DimPh htmp = w0_ph*(iph-1) call sp_insert_element(spH0_ph,htmp,iph,iph) enddo + if(A_ph/=0.d0)then + do iph=1,DimPh + if(iph < DimPh) then !bdg = sum_n |n+1> sqrt(n+1) 1) then !bdg = sum_n |n+1> sqrt(n+1) and with X=(b+bdg)/sqrt(2) + if(iph1 .AND. iph==1) then val = 1 @@ -594,7 +607,7 @@ subroutine write_legend() write(unit,"(A1,90(A10,6X))") "# ",((reg(txtfy(iorb+(jorb-1)*Norb))//"phisc_"//reg(txtfy(iorb)),iorb=1,Norb),jorb=1,Norb) write(unit,"(A1,90(A10,6X))") "# *****" write(unit,"(A1,90(A10,6X))") "# imp_last.ed" - write(unit,"(A1,90(A10,6X))") "# ", "1s2tot", "2egs", "3nph", "4w_ph" + write(unit,"(A1,90(A10,6X))") "# ", "1s2tot", "2egs", "3nph", "4w_ph", "5X_ph", "6X2_ph" close(unit) unit = free_unit() @@ -690,7 +703,7 @@ subroutine write_observables() ! unit = free_unit() open(unit,file="imp_all.ed",position='append') - write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph + write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph, X_ph, X2_ph close(unit) ! !LAST OBSERVABLES @@ -746,7 +759,7 @@ subroutine write_observables() ! unit = free_unit() open(unit,file="imp_last.ed") - write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph + write(unit,"(90(F15.9,1X))") s2tot, egs, dens_ph, w_ph, X_ph, X2_ph close(unit) ! !PARAMETERS diff --git a/src/ED_SUPERC/stored/H_ph.f90 b/src/ED_SUPERC/stored/H_ph.f90 index 9fdfff6..5fd897b 100644 --- a/src/ED_SUPERC/stored/H_ph.f90 +++ b/src/ED_SUPERC/stored/H_ph.f90 @@ -4,3 +4,16 @@ htmp = w0_ph*(iph-1) call sp_insert_element(spH0_ph,htmp,iph,iph) enddo + if(A_ph/=0.d0)then + do iph=1,DimPh + if(iph < DimPh) then !bdg = sum_n |n+1> sqrt(n+1) 1) then !bdg = sum_n |n+1> sqrt(n+1) Date: Tue, 16 May 2023 15:08:37 +0200 Subject: [PATCH 3/4] Small bug fix, cvec instead of dvec --- src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 b/src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 index fcef7d9..c859480 100644 --- a/src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 +++ b/src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90 @@ -183,12 +183,12 @@ subroutine observables_superc() ! and with X=(b+bdg)/sqrt(2) if(iph1 .AND. iph==1) then From abe144da69457a34bd3b0bba82378d9c41381107 Mon Sep 17 00:00:00 2001 From: SamueleGiuli Date: Fri, 26 May 2023 16:33:12 +0200 Subject: [PATCH 4/4] Small Update Add subroutine ed_update_input(name,vals) with input name (a string of char) and an vals (rank-1 array of real(8)) This subroutine can change the values of the internal variables spin_field_x/y/z, exc_field and pair_field, passing their name. This allows to mix standard DMFT with standard MF updating at each loop these internal variables --- src/DMFT_ED.f90 | 1 + src/ED_INPUT_VARS.f90 | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/src/DMFT_ED.f90 b/src/DMFT_ED.f90 index 692f83f..0212781 100644 --- a/src/DMFT_ED.f90 +++ b/src/DMFT_ED.f90 @@ -1,6 +1,7 @@ MODULE DMFT_ED USE ED_INPUT_VARS , only: & ed_read_input , & + ed_update_input,& Norb , & Nspin , & Nbath , & diff --git a/src/ED_INPUT_VARS.f90 b/src/ED_INPUT_VARS.f90 index 1d65576..22efbcd 100644 --- a/src/ED_INPUT_VARS.f90 +++ b/src/ED_INPUT_VARS.f90 @@ -330,6 +330,31 @@ subroutine ed_read_input(INPUTunit) call substring_delete(Hfile,".ed") end subroutine ed_read_input + subroutine ed_update_input(name,vals) + character(len=*) :: name + real(8),dimension(:) :: vals + select case (name) + case default + stop "WRONG NAME ON ED_UPDATE_INPUT" + case ("EXC_FIELD") + if(size(vals)/=4)stop "WRONG SIZE IN ED_UPDATE_EXC_FIELD" + exc_field=vals + case ("PAIR_FIELD") + if(size(vals)/=Norb)stop "WRONG SIZE IN ED_UPDATE_PAIR_FIELD" + pair_field=vals + case ("SPIN_FIELD_X") + if(size(vals)/=Norb)stop "WRONG SIZE IN ED_UPDATE_SPIN_FIELD_X" + spin_field_x=vals + case ("SPIN_FIELD_Y") + if(size(vals)/=Norb)stop "WRONG SIZE IN ED_UPDATE_SPIN_FIELD_Y" + spin_field_y=vals + case ("SPIN_FIELD_Z") + if(size(vals)/=Norb)stop "WRONG SIZE IN ED_UPDATE_SPIN_FIELD_Z" + spin_field_z=vals + end select + + end subroutine ed_update_input +