quantum-espresso/GWW/gww/gww_fit.f90

202 lines
5.8 KiB
Fortran

!
! Copyright (C) 2001-2013 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
PROGRAM gww_fit
USE input_gw, ONLY : input_options, read_input_gww
USE io_global, ONLY : stdout, ionode
USE self_energy_storage
USE expansion
USE energies_gww
USE start_end
USE mp_world, ONLY : mpime
USE para_gww
USE times_gw
USE w_divergence
USE constants, ONLY : RYTOEV
USE io_files, ONLY : prefix,tmp_dir
implicit none
INTEGER, EXTERNAL :: find_free_unit
TYPE(input_options) :: options
TYPE(self_storage) :: ss
TYPE(self_expansion) :: se
TYPE(self_on_real) :: sr
TYPE(quasi_particles) :: qp
TYPE(times_freqs) :: tf
TYPE(gv_time) :: gt
INTEGER :: ii, jj, iun1,iun2,iun,idumm, it
REAL(kind=DP) :: rdumm1,rdumm2,rdumm3,ss_r,ss_i,ddumm
CHARACTER(5) :: nfile
COMPLEX(kind=DP) :: zz, sigmac,dsigmac
REAL(kind=DP) :: offset
INTEGER :: i_homo
INTEGER :: ifil!DEBUG
INTEGER, ALLOCATABLE :: order(:)
!setup MPI environment
call startup
write(stdout,*) 'PROGRAM GWW-FIT: Version 0.2'
!initialize arrays
call initialize_quasi_particle(qp)
call initialize_self_storage(ss)
call initialize_self_on_real(sr)
! read in input structure
call read_input_gww(options)
FLUSH(stdout)
call setup_para_gww(options%n, options%max_i, options%i_min, options%i_max)
! setup time/frequency grid if required
call setup_timefreq(tf,options)
FLUSH(stdout)
! call test_fft(tf)
!read in LDA and HF energies from bands.dat
qp%max_i=options%max_i
qp%nspin=options%nspin
allocate(qp%ene_dft_ks(options%max_i,qp%nspin))
allocate(qp%ene_hf(options%max_i,qp%nspin))
allocate(qp%ene_gw(options%max_i,qp%nspin))
allocate(qp%ene_gw_pert(options%max_i,qp%nspin))
allocate(order(options%max_i))
open( unit=iun, file=trim(tmp_dir)//trim(prefix)//'-'//'bands.dat', status='old',form='formatted')
read(iun,*) i_homo
do ii=1,options%max_i
read(iun,*) idumm,qp%ene_dft_ks(ii,1),rdumm1,rdumm2,qp%ene_hf(ii,1)
enddo
close(iun)
qp%ene_dft_ks(:,1)=qp%ene_dft_ks(:,1)/RYTOEV
qp%ene_hf(:,1)=qp%ene_hf(:,1)/RYTOEV
!setup ss and read in date from fit graphs
ss%ontime=.false.
ss%whole_s=.false.
ss%n=options%n
ss%n_grid_fit=options%n_grid_fit
ss%max_i=options%max_i
ss%i_min=options%i_min
ss%i_max=options%i_max
ss%tau=options%tau
ss%nspin=1
if(options%l_order) then
open( unit=iun, file=trim(tmp_dir)//trim(prefix)//'-'//'order.dat', status='old',form='formatted')
do ii=1,options%max_i
read(iun,*) idumm,order(ii),ddumm
enddo
close(iun)
else
do ii=1,options%max_i
order(ii)=ii
enddo
endif
allocate(ss%diag_freq_fit(ss%max_i,2*ss%n_grid_fit+1,1))
do ii=ss%i_min,ss%i_max
write(nfile,'(5i1)') &
& order(ii)/10000,mod(order(ii),10000)/1000,mod(order(ii),1000)/100,mod(order(ii),100)/10,mod(order(ii),10)
!ifil=ii-4
!write(nfile,'(5i1)') &
! & ifil/10000,mod(ifil,10000)/1000,mod(ifil,1000)/100,mod(ifil,100)/10,mod(ifil,10)
iun1 = find_free_unit()
open( unit=iun1, file=trim(tmp_dir)//trim(prefix)//'-'//'re_on_im'// nfile, status='old',form='formatted')
open( unit=iun2, file=trim(tmp_dir)//trim(prefix)//'-'//'im_on_im'// nfile, status='old',form='formatted')
do jj=-ss%n_grid_fit,ss%n_grid_fit
read(iun1,*) rdumm1,rdumm2,ss_r,rdumm3
read(iun2,*) rdumm1,rdumm2,ss_i,rdumm3
ss%diag_freq_fit(ii,jj+ss%n_grid_fit+1,1)=cmplx(ss_r,ss_i)
enddo
close(iun1)
close(iun2)
enddo
if(options%n_real_axis>=0) then
call create_self_on_real(options, sr)
endif
call create_self_energy_fit( tf, se, ss, options,sr,.false.)
call free_memory_self_storage(ss)
call free_memory_self_on_real(sr)
offset=-(qp%ene_dft_ks(i_homo+1,1)+qp%ene_dft_ks(i_homo,1))/2.d0
!offset=-(6.07889+6.75583)/2.d0/13.606d0!DEBUG
do ii=options%i_min,options%i_max
call value_on_frequency(se,ii,qp%ene_dft_ks(ii,1)+offset,sigmac,1)
call derivative_on_frequency(se,ii,qp%ene_dft_ks(ii,1)+offset,dsigmac,1)
write(stdout,*) 'value, zeta:',ii,sigmac,dsigmac,offset
zz=(1.d0,0.d0)-dsigmac
qp%ene_gw(ii,1)=qp%ene_dft_ks(ii,1)+offset +&
(sigmac+qp%ene_hf(ii,1)-qp%ene_dft_ks(ii,1))/zz
write(stdout,*) 'GW-PERT energy', ii,real(qp%ene_gw(ii,1)-offset)*RYTOEV
qp%ene_gw_pert(ii,1)=qp%ene_gw(ii,1)-offset
!self-consistency loop
do it=1,10
call value_on_frequency_complex(se,ii,qp%ene_gw(ii,1),sigmac,1)
write(stdout,*) 'Iteration energy',it,sigmac
qp%ene_gw(ii,1)=qp%ene_dft_ks(ii,1)+offset+sigmac+qp%ene_hf(ii,1)-qp%ene_dft_ks(ii,1)
enddo
qp%ene_gw(ii,1)= qp%ene_gw(ii,1)-offset
enddo
!call create_quasi_particles(options,qp,se)
call free_memory_self_expansion(se)
if(ionode) then
write(stdout,*) 'QUASI-PARTICLES ENERGIES IN Ev:'
do ii=options%i_min,options%i_max
write(stdout,'(''State:'',i5,''LDA :'',f10.5,'' GW-PERT :'',f10.5,'' GW :'',f10.5, &
& '' HF-pert :'',f10.5)') &
& ii,qp%ene_dft_ks(ii,1)*RYTOEV, real(qp%ene_gw_pert(ii,1))*RYTOEV, &
& real(qp%ene_gw(ii,1))*RYTOEV,qp%ene_hf(ii,1)*RYTOEV
enddo
write(stdout,*) 'IMAGINARY ENERGIES IN Ev:'
do ii=options%i_min,options%i_max
write(stdout,'(''State:'',i5,'' GW (Im) :'',f10.5)') ii,aimag(qp%ene_gw(ii,1))*RYTOEV
enddo
endif
deallocate(order)
!stops MPI
call free_memory_times_freqs(tf)
call free_memory_para_gww
call stop_run
stop
END PROGRAM gww_fit