!
! License-Identifier: GPL
!
! Copyright (C) 2016 The Yambo Team
!
! Authors (see AUTHORS file for details): AM DV
!
! headers
!
#include<y_memory.h>
!
subroutine plot_gnuplot(eval_only)
!
use pars, ONLY:SP,schlen
use units, ONLY:BO2ANG
use com, ONLY:msg
use YPPm, ONLY:nr,v2plot,r_fixed,v2plot2D,l_norm_to_one,WF_multiplier,plot_dim,&
& plot_title,output_string,l_wavefunction,l_current,use_xcrysden,&
& plot_is_1D,plot_is_2D,l_free_hole,l_average_eh
use LIVE_t, ONLY:live_timing
use y_memory_alloc
!
implicit none
logical :: eval_only
!
! Work Space...
!
integer :: i1,i2,i3,ir,nr2plot(2),dir2plot(2)
character(schlen) :: titles(4)
real(SP) :: v_max,r_plot_center(3)
!
if (plot_dim==3) return
!
if (len_trim(plot_title)==0) then
call live_timing('3D Merge',nr(1))
else
call live_timing('3D Merge of '//trim(plot_title),nr(1))
endif
!
if (.not.eval_only) then
!
titles=" "
!
titles(:3)=(/'||a1 [A]','||a2 [A]','||a3 [A]'/)
if (l_wavefunction) titles(4) = '|wf|^2 [1 at max]'
if (l_current) titles(4) = 'J [nAmpere]'
!
endif
!
! SET CENTER OF THE PLOT
!
r_plot_center=0.
if (.not.l_free_hole.and..not.l_average_eh) r_plot_center=r_fixed
!
! DIMENSIONs
!
select case(plot_dim)
!
case(1)
!
do i1=1,3
if (plot_is_1D(i1)) nr2plot=nr(i1)
if (plot_is_1D(i1)) dir2plot=i1
if (plot_is_1D(i1)) call msg(output_string,"#",(/titles(1),titles(4)/),INDENT=0,USE_TABS=.true.)
enddo
call msg(output_string,'#')
YAMBO_ALLOC(v2plot2D,(nr2plot(1),1))
v2plot2D=0._SP
ir=0
do i1 = 0,nr(1)-1
do i2 = 0,nr(2)-1
do i3 = 0,nr(3)-1
ir = 1 + i1 + i2*nr(1) + i3*nr(1)*nr(2)
if (plot_is_1D(1)) v2plot2D(i1+1,1)=v2plot2D(i1+1,1)+v2plot(ir)
if (plot_is_1D(2)) v2plot2D(i2+1,1)=v2plot2D(i2+1,1)+v2plot(ir)
if (plot_is_1D(3)) v2plot2D(i3+1,1)=v2plot2D(i3+1,1)+v2plot(ir)
enddo
enddo
call live_timing(steps=1)
enddo
v_max=maxval(v2plot2D(:,1))
!
if (l_norm_to_one) v2plot2D=v2plot2D/v_max
v2plot2D=v2plot2D*WF_multiplier
!
do i1=0,nr2plot(1)-1
call msg(output_string,'',(/projection_(i1,dir2plot(1)),&
& v2plot2D(i1+1,1)/),INDENT=-2,USE_TABS=.true.)
enddo
!
YAMBO_FREE(v2plot2D)
!
case(2)
!
if (.not.use_xcrysden) then
if (plot_is_2D(1)) nr2plot=(/nr(1),nr(2)/)
if (plot_is_2D(2)) nr2plot=(/nr(1),nr(3)/)
if (plot_is_2D(3)) nr2plot=(/nr(2),nr(3)/)
else
if (plot_is_2D(1)) nr2plot=(/nr(1)+1,nr(2)+1/)
if (plot_is_2D(2)) nr2plot=(/nr(1)+1,nr(3)+1/)
if (plot_is_2D(3)) nr2plot=(/nr(2)+1,nr(3)+1/)
endif
if (plot_is_2D(1)) dir2plot=(/1,2/)
if (plot_is_2D(2)) dir2plot=(/1,3/)
if (plot_is_2D(3)) dir2plot=(/2,3/)
YAMBO_ALLOC(v2plot2D,(nr2plot(1),nr2plot(2)))
v2plot2D=0.
ir=0
do i1 = 0,nr(1)-1
do i2 = 0,nr(2)-1
do i3 = 0,nr(3)-1
ir = 1 + i1 + i2*nr(1) + i3*nr(1)*nr(2)
if (plot_is_2D(1)) v2plot2D(i1+1,i2+1)=v2plot2D(i1+1,i2+1)+v2plot(ir)
if (plot_is_2D(2)) v2plot2D(i1+1,i3+1)=v2plot2D(i1+1,i3+1)+v2plot(ir)
if (plot_is_2D(3)) v2plot2D(i2+1,i3+1)=v2plot2D(i2+1,i3+1)+v2plot(ir)
enddo
enddo
call live_timing(steps=1)
enddo
if(use_xcrysden) then
v2plot2D(nr2plot(1),:nr2plot(2)-1)=v2plot2D(1,:nr2plot(2)-1)
v2plot2D(:nr2plot(1)-1,nr2plot(2))=v2plot2D(:nr2plot(1)-1,1)
v2plot2D(nr2plot(1),nr2plot(2))=v2plot2D(1,1)
endif
!
v_max=maxval(v2plot2D(:,:))
if (l_norm_to_one) v2plot2D=v2plot2D/v_max
v2plot2D=v2plot2D*WF_multiplier
!
if (.not.eval_only) then
if (plot_is_2D(1)) call msg(output_string,"#",&
& (/titles(1),titles(2),titles(4)/),INDENT=0,USE_TABS=.true.)
if (plot_is_2D(2)) call msg(output_string,"#",&
& (/titles(1),titles(3),titles(4)/),INDENT=0,USE_TABS=.true.)
if (plot_is_2D(3)) call msg(output_string,"#",&
& (/titles(2),titles(3),titles(4)/),INDENT=0,USE_TABS=.true.)
call msg(output_string,"#")
do i1=0,nr2plot(1)-1
do i2=0,nr2plot(2)-1
call msg(output_string,'',(/projection_(i1,dir2plot(1)),&
& projection_(i2,dir2plot(2)),v2plot2D(i1+1,i2+1)/),&
& INDENT=-2,USE_TABS=.true.)
enddo
enddo
endif
!
if (.not.use_xcrysden) then
!
YAMBO_FREE(v2plot2D)
!
endif
!
end select
!
call live_timing()
!
contains
!
real(SP) function projection_(I,dir)
use D_lattice, ONLY:a
use FFT_m, ONLY:fft_dim
use vec_operate, ONLY:v_norm
integer :: I,dir
real(SP):: rv(3)
!
rv(:)=I*a(dir,:)/fft_dim(dir)-r_plot_center(:)
projection_=dot_product(rv,a(dir,:))/v_norm(a(dir,:))*BO2ANG
!
end function
!
end subroutine
