Merge branch 'master' of github.com:aamaricci/EDIpack2.0

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 , &

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))
@@ -328,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
+
 
 
 

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)

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
  !

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
+
+ !<X> and <X^2> with X=(b+bdg)/sqrt(2)
+ if(iph<DimPh)then
+ j= i_el + (iph)*sectorI%DimEl
+ X_ph = X_ph + sqrt(2.d0*dble(iph))*(state_dvec(i)*state_dvec(j))*peso
+ end if
+ X2_ph = X2_ph + 0.5d0*(1+2*(iph-1))*gs_weight
+ if(iph<DimPh-1)then
+ j= i_el + (iph+1)*sectorI%DimEl
+ X2_ph = X2_ph + sqrt(dble((iph)*(iph+1)))*(state_dvec(i)*state_dvec(j))*peso
+ end if
  !
  !compute the lattice probability distribution function
  if(Dimph>1 .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()

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

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)<n|
+ htmp = sqrt(0.5d0)*A_ph*sqrt(dble(iph))
+ call sp_insert_element(spH0_ph,htmp,iph+1,iph)
+ endif
+ if(iph > 1) then !bdg = sum_n |n+1> sqrt(n+1)<n|
+ htmp = sqrt(0.5d0)*A_ph*sqrt(dble(iph-1))
+ call sp_insert_element(spH0_ph,htmp,iph-1,iph)
+ endif
+ enddo
+ end if

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

src/ED_SUPERC/ED_OBSERVABLES_SUPERC.f90

@@ -35,6 +35,7 @@ MODULE ED_OBSERVABLES_SUPERC
  real(8),dimension(:),allocatable :: pdf_ph
  real(8),dimension(:,:),allocatable :: pdf_part
  real(8) :: dens_ph
+ real(8) :: X_ph, X2_ph
  real(8) :: w_ph
  !
  integer :: iorb,jorb,istate
@@ -103,6 +104,8 @@ subroutine observables_superc()
  s2tot = 0.d0
  Prob = 0.d0
  prob_ph = 0.d0
+ X_ph = 0.d0
+ X2_ph = 0.d0
  dens_ph = 0.d0
  pdf_ph = 0.d0
  pdf_part= 0.d0
@@ -177,6 +180,16 @@ subroutine observables_superc()
  prob_ph(iph) = prob_ph(iph) + gs_weight
  dens_ph = dens_ph + (iph-1)*gs_weight
  !
+ !<X> and <X^2> with X=(b+bdg)/sqrt(2)
+ if(iph<DimPh)then
+ j= i_el + (iph)*sectorI%DimEl
+ X_ph = X_ph + sqrt(2.d0*dble(iph))*real(state_cvec(i)*conjg(state_cvec(j)))*peso
+ end if
+ X2_ph = X2_ph + 0.5d0*(1+2*(iph-1))*gs_weight
+ if(iph<DimPh-1)then
+ j= i_el + (iph+1)*sectorI%DimEl
+ X2_ph = X2_ph + sqrt(dble((iph)*(iph+1)))*real(state_cvec(i)*conjg(state_cvec(j)))*peso
+ end if
  !compute the lattice probability distribution function
  if(Dimph>1 .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

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)<n|
+ htmp = sqrt(0.5d0)*A_ph*sqrt(dble(iph))
+ call sp_insert_element(spH0_ph,htmp,iph+1,iph)
+ endif
+ if(iph > 1) then !bdg = sum_n |n+1> sqrt(n+1)<n|
+ htmp = sqrt(0.5d0)*A_ph*sqrt(dble(iph-1))
+ call sp_insert_element(spH0_ph,htmp,iph-1,iph)
+ endif
+ enddo
+ end if
+