diff --git a/src/bse/K_driver_init.F b/src/bse/K_driver_init.F index 0ace8836ab..1e73150fb1 100644 --- a/src/bse/K_driver_init.F +++ b/src/bse/K_driver_init.F @@ -28,7 +28,7 @@ subroutine K_driver_init(what,iq,Ken,Xk) use drivers, ONLY:l_rt_carriers_in_use use RT_control, ONLY:NEQ_Kernel,EQ_Transitions,EQ_NoOcc,NEQ_Residuals,RT_BSE_Occ_Mode #endif -#if defined _CUDA +#if defined _GPU use TDDFT, ONLY:FXC_mode #endif ! @@ -173,7 +173,7 @@ subroutine K_driver_init(what,iq,Ken,Xk) call parser('FxcLibxc',l_Fxc_Libxc) l_Fxc_from_Vxc=(n_spin==1).and..not.l_Fxc_Libxc endif -#if defined _CUDA +#if defined _GPU if ( BS_K_is_alda.and. index(FXC_mode,"G-")>0 ) then call warning(" Tddft with G-integrals is not GPU ported. Fallback to R-integrals") FXC_mode="R-def" diff --git a/src/bse/K_exchange_collisions.F b/src/bse/K_exchange_collisions.F index ef1bc31cbc..1e1de3ba75 100644 --- a/src/bse/K_exchange_collisions.F +++ b/src/bse/K_exchange_collisions.F @@ -82,15 +82,15 @@ subroutine K_exchange_collisions(iq,Xk,i_T_grp,NG,l_bs_exch_wf_in_loop) ! l_load_WFs= l_bs_exch_wf_in_loop .and. (NK(1)/=min(i_k,i_p).or.NK(2)/=max(i_k,i_p)) if (l_load_WFs) then -#if defined(__NOTNOW) && ! defined(_CUDA) - !$omp critical +#if defined(__NOTNOW) + !DEV_OMP critical #endif if (NK(2)/=-1) call WF_free(WF,keep_fft=.true.,keep_states_to_load=.true.) NK=(/min(i_k,i_p),max(i_k,i_p)/) call WF_load(WF,NG(1),NG(2),BS_bands,NK,k_extrema_only=.true.,quiet=.true.,& & space='R',title="Kernel exch",keep_states_to_load=.true.) -#if defined(__NOTNOW) && ! defined(_CUDA) - !$omp end critical +#if defined(__NOTNOW) + !DEV_OMP end critical #endif endif ! diff --git a/src/bse/K_kernel.F b/src/bse/K_kernel.F index 0a66333988..66a782e9c2 100644 --- a/src/bse/K_kernel.F +++ b/src/bse/K_kernel.F @@ -567,7 +567,6 @@ subroutine K_kernel(iq,Ken,Xk,q,X,Xw,W_bss) ! call OPENMP_update(master_thread) ! - ! ! Allocate tddft_wf variables neeeded within the OMP loop if (l_tddft_rsum.and.iHxc==2) then YAMBO_ALLOC(tddft_wf%WF_symm1,(fft_size,n_spinor)) @@ -709,11 +708,14 @@ subroutine K_kernel(iq,Ken,Xk,q,X,Xw,W_bss) ! if (iHxc==3) then ! - if (iq_W_bz/=iq_W_bz_mq.or.iq_W/=iq_W_mq.or.iq_W_s/=iq_W_s_mq) call error("Wrong transferred momentum") + if (iq_W_bz/=iq_W_bz_mq.or.iq_W/=iq_W_mq.or.iq_W_s/=iq_W_s_mq) & +& call error("Wrong transferred momentum") ! - !if ( G_m_G(ig_W,ig_W_mq) /= G_m_G(ig_kmq,ig_pmq) ) call error("Wrong gvector shifts") + !if ( G_m_G(ig_W,ig_W_mq) /= G_m_G(ig_kmq,ig_pmq) ) & +!& call error("Wrong gvector shifts") ! - if ( (.not.BS_W_is_diagonal) .and. iq_W_s>nsym/(i_time_rev+1) .and. i_space_inv == 0 ) iq_W=q%nibz+iq_W + if ( (.not.BS_W_is_diagonal) .and. iq_W_s>nsym/(i_time_rev+1) & +& .and. i_space_inv == 0 ) iq_W=q%nibz+iq_W ! endif ! diff --git a/src/coulomb/rim.F b/src/coulomb/rim.F index b958b4af11..faaa20cedb 100644 --- a/src/coulomb/rim.F +++ b/src/coulomb/rim.F @@ -3,7 +3,9 @@ ! ! Copyright (C) 2006 The Yambo Team ! -! Authors (see AUTHORS file for details): AM +! Authors (see AUTHORS file for details): AM AF NS GS +! +#include ! subroutine rim(mode,Xw) ! @@ -17,10 +19,12 @@ subroutine rim(mode,Xw) use R_lattice, ONLY:RL_vol,k_grid_uc_vol,k_grid_b,nqbz,& & nqibz,RIM_epsm1,RIM_is_diagonal,RIM_RL_vol,& & RIM_n_rand_pts,RIM_ng,RIM_W_ng,RIM_qpg,RIM_W_is_diagonal,& -& RIM_id_epsm1_reference,RIM_anisotropy,RIM_W,& -& cut_is_slab,idir +& RIM_id_epsm1_reference,RIM_anisotropy,RIM_W,RIM_W_d,aux_RIM_W,& ! aux_RIM_W_d,& +& cut_is_slab,idir,f_coeff,f_coeff_d use timing_m, ONLY:timing use frequency, ONLY:w_samp + use gpu_m, ONLY:have_gpu + use devxlib, ONLY:devxlib_memcpy_h2d ! #include ! @@ -29,9 +33,9 @@ subroutine rim(mode,Xw) ! ! Work Space ! - type(PP_indexes) ::px - integer :: iq - real(SP) :: em1_anis(3),G_radii,G_circ + type(PP_indexes) :: px + integer :: iq + real(SP) :: em1_anis(3),G_radii,G_circ ! ! Random generator ! @@ -39,9 +43,10 @@ subroutine rim(mode,Xw) integer :: N_out,N_in,N_out_G,rep_factor,inn1,inn2,inn3,ic real(SP) :: v1(3),v2(3),v1_norm(2) real(SP), allocatable :: qr(:,:) + real(SP), allocatable DEV_ATTR :: qr_d(:,:) integer :: iseed(8) real(DP), external :: dlaran - ! + ! if (mode == "x") call timing('RIM',OPR='start') if (mode == "c") call timing('RIM-W',OPR='start') ! @@ -53,6 +58,9 @@ subroutine rim(mode,Xw) ! the bare part is embodied in the exchange. ! YAMBO_ALLOC(qr,(3,RIM_n_rand_pts)) + if (have_gpu) then + YAMBO_ALLOC_GPU(DEV_VAR(qr),(3,RIM_n_rand_pts)) + endif ! em1_anis=RIM_epsm1(:)-1. ! @@ -146,6 +154,10 @@ subroutine rim(mode,Xw) !Only the 2D-BZ is sampled for 2D systems if (cut_is_slab) qr(idir(1),:) = 0._SP ! + ! init on GPU mem + ! + if (have_gpu) call devxlib_memcpy_h2d(DEV_VAR(qr),qr) + ! call live_timing() call msg('r','Points outside the sBZ ',N_out) ! @@ -166,7 +178,16 @@ subroutine rim(mode,Xw) else if (mode == "c") then YAMBO_ALLOC(RIM_W,(Xw%n_freqs,nqibz,RIM_W_ng,RIM_W_ng)) RIM_W=0._SP - end if + ! + ! AF: need to cp f_coeff on GPU mem before calling rim_integrate_w + ! + YAMBO_ALLOC(aux_RIM_W,(Xw%n_freqs,RIM_W_ng,RIM_W_ng)) +#ifdef _GPU + !AMBO_ALLOC_GPU_SOURCE(DEV_VAR(aux_RIM_W),aux_RIM_W) + YAMBO_ALLOC_GPU_SOURCE(DEV_VAR(RIM_W),RIM_W) + YAMBO_ALLOC_GPU_SOURCE(DEV_VAR(f_coeff),f_coeff) +#endif + endif ! call PARALLEL_index(px,(/nqibz/)) call live_timing('Momenta loop',px%n_of_elements(myid+1)) @@ -176,8 +197,10 @@ subroutine rim(mode,Xw) ! if (mode == "x") then call rim_integrate_v(iq,qr,em1_anis,N_out,N_out_G,G_radii,G_circ) - else - call rim_integrate_w(iq,qr,N_out,em1_anis,Xw) + else if (mode == "c") then + ! + call rim_integrate_w(iq,qr,qr_d,N_out,em1_anis,Xw) + ! endif call live_timing(steps=1) ! @@ -210,7 +233,16 @@ subroutine rim(mode,Xw) ! CLEAN ! call PP_indexes_reset(px) + ! + if (mode == "c") then + YAMBO_FREE(aux_RIM_W) +#ifdef _GPU + YAMBO_FREE_GPU(DEV_VAR(f_coeff)) + YAMBO_FREE_GPU(DEV_VAR(qr)) + !AMBO_FREE_GPU(DEV_VAR(aux_RIM_W)) +#endif YAMBO_FREE(qr) + end if ! if (mode == "x") call timing('RIM',OPR='stop') if (mode == "c") call timing('RIM-W',OPR='stop') diff --git a/src/coulomb/rim_integrate_w.F b/src/coulomb/rim_integrate_w.F index 5f0cfcb7c5..8981a21305 100644 --- a/src/coulomb/rim_integrate_w.F +++ b/src/coulomb/rim_integrate_w.F @@ -1,156 +1,986 @@ ! ! License-Identifier: GPL +! Copyright (C) 2000-2022 the YAMBO team +! http://www.yambo-code.org ! ! Copyright (C) 2006 The Yambo Team ! -! Authors (see AUTHORS file for details): AM +#include ! -subroutine rim_integrate_w(iq,qr,N_out,em1_anis,Xw) +! Authors (see AUTHORS file for details): AM AG GS AF NS +! +! AF: This file needs cleanup +! - proper indentation +! - explicit handling of precision (e.g. 2.0_SP instead of 2.0) +! - too long, splitting ?? +! +subroutine rim_integrate_w(iq,qr,qr_d,N_out,em1_anis,Xw) ! - use pars, ONLY:SP,DP,pi - use vec_operate, ONLY:iku_v_norm - use R_lattice, ONLY:g_vec,RIM_n_rand_pts,k_grid_uc_vol,q0_def_norm,& -& RIM_W_ng,q_pt,b,RIM_W_is_diagonal,& -& RIM_W,f_coeff,cut_is_slab,idir,RIM_id_epsm1_reference + use pars, ONLY:SP,DP,pi,zero_dfl + use com, ONLY:msg + use vec_operate, ONLY:c2a,iku_v_norm + use R_lattice, ONLY:g_vec,RIM_n_rand_pts,k_grid_b,k_grid_uc_vol,q0_def_norm,& +& RIM_W_ng,q_pt,b,RIM_W_is_diagonal,idir,RIM_W,RIM_W_d,& +& aux_RIM_W,f_coeff,DEV_VAR(f_coeff),cut_is_slab,RIM_id_epsm1_reference,rimw_type use D_lattice, ONLY:alat,a use frequency, ONLY:w_samp + use gpu_m, ONLY:have_gpu + use devxlib, ONLY:devxlib_memcpy_d2h,devxlib_memcpy_h2d + ! +#include ! - implicit none - integer :: iq,N_out - real(SP) :: qr(3,RIM_n_rand_pts),em1_anis(3) + integer :: iq,N_out,i_dir + real(SP) :: em1_anis(3),qtmp(3),qr_rlu(3) + real(SP) :: qr(3,RIM_n_rand_pts) + real(SP) DEV_ATTR :: qr_d(3,RIM_n_rand_pts) + real(SP) :: v1(3),q_grid_b_rlu(3,3),q_grid_b_iku(3,3) type(w_samp) :: Xw ! ! Work Space ! - integer :: i1,i1min,i2,i2max,i3,iw - real(DP) :: vslab(RIM_n_rand_pts,RIM_W_ng) + integer :: i1,i1min,i2,i2min,i2max,i3,iw,j + real(DP), allocatable :: vslab(:,:) + real(DP), allocatable DEV_ATTR :: vslab_d(:,:) real(DP) :: r1,rfac complex(DP):: RIM_acc,RIM_acc_anis,func real(SP) :: slab_vz1,slab_vplane1,lcut,pre_factor,anis_fact,vslab2 ! + ! define pointers to enable CUF kernels + ! when compiling using CUDA-Fortran + ! + complex(DP), pointer DEV_ATTR :: f_coeff_p(:,:,:,:,:) + ! + ! AF: f_coeff_p needs to be present on GPU memory + f_coeff_p => DEV_VAR(f_coeff) + ! + call msg("r","start rimw integration") + call msg("r", trim(rimw_type)) + ! rfac=8._SP*k_grid_uc_vol/real(N_out)/(2._SP*pi)**4 ! if (cut_is_slab) then - lcut=alat(idir(1))/2._SP - else - call error('RIM-W without cutoff slab has not been implemented') + lcut=alat(idir(1))/2._SP + ! else + ! call error('RIM-W without cutoff slab has not been implemented') end if ! - !Evaluate vslab + !Load q_grid_b + do j=1,3 + v1 = k_grid_b(j,:) + call c2a(v_in=v1,mode='kc2a') + ! + if ((abs(v1(2))>> do i3=1,RIM_n_rand_pts ! - func = f_coeff(1,i1,i2,iq,iw)+qr(1,i3)*(f_coeff(2,i1,i2,iq,iw)+qr(1,i3)*f_coeff(4,i1,i2,iq,iw))& -& +qr(2,i3)*(f_coeff(3,i1,i2,iq,iw)+qr(2,i3)*f_coeff(6,i1,i2,iq,iw)& -& +2._DP*qr(1,i3)*f_coeff(5,i1,i2,iq,iw)) - ! Accumulate W - RIM_acc = RIM_acc + rfac*vslab(i3,i1)**2*func*vslab(i3,i2)**2/(1._DP-vslab(i3,i1)*func*vslab(i3,i2)) + func = f_coeff_p(1,i1,i2,iq,iw) & +& + DEV_VAR(qr)(1,i3)*(f_coeff_p(2,i1,i2,iq,iw)+DEV_VAR(qr)(1,i3)*f_coeff_p(5,i1,i2,iq,iw) & +& + 2._DP*DEV_VAR(qr)(2,i3)*f_coeff_p(6,i1,i2,iq,iw)) & +& + DEV_VAR(qr)(2,i3)*(f_coeff_p(3,i1,i2,iq,iw)+DEV_VAR(qr)(2,i3)*f_coeff_p(8,i1,i2,iq,iw) & +& + 2._DP*DEV_VAR(qr)(3,i3)*f_coeff_p(9,i1,i2,iq,iw)) & +& + DEV_VAR(qr)(3,i3)*(f_coeff_p(4,i1,i2,iq,iw)+DEV_VAR(qr)(3,i3)*f_coeff_p(10,i1,i2,iq,iw) & +& + 2._DP*DEV_VAR(qr)(1,i3)*f_coeff_p(7,i1,i2,iq,iw)) & +& +DEV_VAR(qr)(1,i3)**2*(DEV_VAR(qr)(1,i3)*f_coeff_p(11,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(2,i3)* & +& f_coeff_p(12,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(3,i3)*f_coeff_p(13,i1,i2,iq,iw)) & +& +DEV_VAR(qr)(2,i3)**2*(DEV_VAR(qr)(2,i3)*f_coeff_p(15,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(3,i3)* & +& f_coeff_p(16,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(1,i3)*f_coeff_p(14,i1,i2,iq,iw)) & +& +DEV_VAR(qr)(3,i3)**2*(DEV_VAR(qr)(3,i3)*f_coeff_p(19,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(1,i3)* & +& f_coeff_p(17,i1,i2,iq,iw)+3_DP*DEV_VAR(qr)(2,i3)*f_coeff_p(18,i1,i2,iq,iw)) & +& +6*DEV_VAR(qr)(1,i3)*DEV_VAR(qr)(2,i3)*DEV_VAR(qr)(3,i3)*f_coeff_p(20,i1,i2,iq,iw) + ! + ! Accumulate Wc + ! + RIM_acc = RIM_acc + & +& rfac*DEV_VAR(vslab)(i3,i1)**2*func*DEV_VAR(vslab)(i3,i2)**2 / & +& (1._DP- DEV_VAR(vslab)(i3,i1)*func*DEV_VAR(vslab)(i3,i2) ) + ! enddo ! - RIM_W(iw,iq,i1,i2)=RIM_acc - RIM_W(iw,iq,i2,i1)=RIM_W(iw,iq,i1,i2) + aux_RIM_W(iw,i1,i2)=cmplx(RIM_acc,KIND=SP) ! enddo enddo + ! + !DEV_OMP end parallel do + ! enddo - !$omp end parallel do + ! + !if (have_gpu) call devxlib_memcpy_d2h(aux_RIM_W,DEV_VAR(aux_RIM_W)) + ! + RIM_W(:,iq,:,:)=aux_RIM_W(:,:,:) + ! +#ifdef _GPU + YAMBO_FREE_GPU(DEV_VAR(vslab)) +#endif + ! if (iq>1) return ! ! head q == 0 ! do iw=1,Xw%n_freqs + ! RIM_acc=0._DP RIM_acc_anis=0._DP ! - do i1=1,RIM_n_rand_pts - ! - slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2+& -& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) - !kxy - r1=iku_v_norm(qr(:,i1)) - !Regularization - if (slab_vplane1 < 1.e-5) then - vslab(i1,1) = sqrt(4._DP*pi*(1.-exp(-q0_def_norm*lcut))/q0_def_norm**2) - RIM_acc = RIM_acc + rfac*f_coeff(1,1,1,1,iw)*(4*pi*lcut)**2 - RIM_acc_anis = RIM_acc_anis + rfac*f_coeff(1,1,1,1,iw)*(4*pi*lcut)**2 & -& *0.5_SP*(em1_anis(idir(2))+em1_anis(idir(3))) - cycle - end if + if (iw==1) then + select case (trim(rimw_type)) + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("semiconductor") + do i1=1,RIM_n_rand_pts + ! + call c2a(v_in=qr(:,i1),v_out=qr_rlu(:),mode='ki2a') + ! + slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) + !kxy + r1=sqrt((qr(1,i1)*2.*pi/alat(1))**2+(qr(2,i1)*2.*pi/alat(2))**2 + & +& (qr(3,i1)*2.*pi/alat(3))**2) + ! + !Regularization + if (r1 < 1.e-5) then + ! + if (cut_is_slab) then + ! + ! System 2D, prefactor of slab truncation + ! + vslab(i1,1) = sqrt(4._DP*pi*(1.-exp(-q0_def_norm*lcut))/q0_def_norm**2) + RIM_acc= RIM_acc + rfac*f_coeff(1,1,1,1,iw)*(4*pi*lcut)**2 + RIM_acc_anis= RIM_acc_anis + rfac*f_coeff(1,1,1,1,iw)*(4.0_DP*pi*lcut)**2 & +& *0.5_DP*(em1_anis(idir(2))+em1_anis(idir(3))) + ! + else + ! + ! System 3D, standard bare coulomb + ! + vslab(i1,1) = sqrt(4._DP*pi/q0_def_norm**2) + RIM_acc= RIM_acc +rfac*f_coeff(1,1,1,1,iw)/(1.0-(4._DP*pi)*f_coeff(1,1,1,1,iw)) & +& *(4._DP*pi)**2/q0_def_norm**2 + RIM_acc_anis= RIM_acc_anis + rfac*f_coeff(1,1,1,1,iw)/(1.0_DP-(4._DP*pi) & + *f_coeff(1,1,1,1,iw))*(4._DP*pi)**2/q0_def_norm**2 & +& *(em1_anis(1)+em1_anis(2)+em1_anis(3))/3.0_SP + ! + endif + ! + cycle + endif + ! + !Evaluate v_slab and interpolation function + ! + if (cut_is_slab) then + ! + vslab2 = 4._DP*pi*(1.-exp(-slab_vplane1*lcut)) + func = f_coeff(1,1,1,1,iw)*exp(sign(1._DP,real(f_coeff(idir(2)+1,1,1,1,iw),DP)) & +& *sqrt((f_coeff(idir(2)+1,1,1,1,iw)*qr_rlu(idir(2)))**2 & +& +(f_coeff(idir(3)+1,1,1,1,iw)*qr_rlu(idir(3)))**2)) + ! + else + ! + vslab2 = 4._DP*pi + func = f_coeff(1,1,1,1,iw)*exp(sign(1._DP,real(f_coeff(2,1,1,1,iw),DP)) & +& *sqrt((f_coeff(2,1,1,1,iw)*qr_rlu(1))**2 + (f_coeff(3,1,1,1,iw) & +& *qr_rlu(2))**2 + (f_coeff(4,1,1,1,iw)*qr_rlu(3))**2)) + ! + end if + ! + !Evaluate W + RIM_acc = RIM_acc + rfac*vslab2*func*vslab2/(1-vslab2*func)/r1**2 + ! + !Anisotropy contribution + anis_fact=dot_product(em1_anis,(2._SP*pi*qr(:,i1)/(alat(:)*r1))**2) + func = func*anis_fact + RIM_acc_anis = RIM_acc_anis + rfac*vslab2*func*vslab2/(r1**2*(1-vslab2*func)) + ! + !Store the square root of vslab + vslab(i1,1) = sqrt(vslab2)/r1 + ! + enddo + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("Dirac") + do i1=1,RIM_n_rand_pts + ! + call c2a(v_in=qr(:,i1),v_out=qr_rlu(:),mode='ki2a') + ! + slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) + !kxy + r1=sqrt((qr(1,i1)*2.*pi/alat(1))**2+(qr(2,i1)*2.*pi/alat(2))**2 + & +& (qr(3,i1)*2.*pi/alat(3))**2) + ! + if (r1 < 1.e-5) then + vslab2 = 4._DP*pi/q0_def_norm**2 + ! + ! System 2D, prefactor of slab truncation + ! + if (cut_is_slab) vslab2 = vslab2*sqrt(1.-exp(-q0_def_norm*lcut)) + ! + func = f_coeff(1,1,1,1,1)*q0_def_norm + RIM_acc = RIM_acc + rfac*vslab(i1,1)**2*func*vslab(i1,1)**2 & +& /(1._DP-vslab(i1,1)**2*func) + cycle + ! + endif + ! + !Evaluate v_slab + vslab2 = 4._DP*pi/r1**2 + ! + ! System 2D, prefactor of slab truncation + ! + if (cut_is_slab) vslab2 = vslab2*(1.-exp(-slab_vplane1*lcut)) + ! + !Evaluate func + ! + func = f_coeff(1,1,1,1,iw)*r1*exp(sign(1._DP,real(f_coeff(idir(2)+1,1,1,1,iw),DP)) & +& *sqrt((f_coeff(idir(2)+1,1,1,1,iw)*qr_rlu(idir(2)))**2 & +& +(f_coeff(idir(3)+1,1,1,1,iw)*qr_rlu(idir(3)))**2)) + ! + RIM_acc = RIM_acc+rfac*vslab2*func*vslab2/(1._DP-vslab2*func) + ! + !Store the square root of vslab + vslab(i1,1) = sqrt(vslab2) + ! + enddo + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("metal") + do i1=1,RIM_n_rand_pts + ! + call c2a(v_in=qr(:,i1),v_out=qr_rlu(:),mode='ki2a') + ! + func = f_coeff(1,1,1,1,1) + f_coeff(2,1,1,1,1)*abs(qr_rlu(1)) + f_coeff(3,1,1,1,1) & +& *abs(qr_rlu(2)) + f_coeff(4,1,1,1,1)*abs(qr_rlu(3)) + f_coeff(5,1,1,1,1) & +& *qr_rlu(1)**2 + f_coeff(6,1,1,1,1)*qr_rlu(2)**2 + f_coeff(7,1,1,1,1)*qr_rlu(3)**2 + ! + slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2+& +& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) + ! + r1=sqrt((qr(1,i1)*2.*pi/alat(1))**2+(qr(2,i1)*2.*pi/alat(2))**2 + & +& (qr(3,i1)*2.*pi/alat(3))**2) + ! + if (r1 < 1.e-5) then + ! + vslab(i1,1) = sqrt(4._DP*pi/q0_def_norm**2) + ! + if (cut_is_slab) vslab(i1,1) = vslab(i1,1)*sqrt(1.-exp(-q0_def_norm*lcut)) + ! + else + ! + vslab(i1,1) =sqrt(4._DP*pi/r1**2) + ! + if (cut_is_slab) vslab(i1,1) = vslab(i1,1)*sqrt(1.-exp(-slab_vplane1*lcut)) + ! + endif + ! + RIM_acc = RIM_acc + rfac*func*vslab(i1,1)**4/(1._DP-func*vslab(i1,1)**2) + ! + enddo ! - !Evaluate v_slab - vslab2=4._DP*pi*(1.-exp(-slab_vplane1*lcut)) + end select ! - !Evaluate interpolation function - func = f_coeff(1,1,1,1,iw)*exp(-sqrt((f_coeff(2,1,1,1,iw)*(qr(2,i1)+a(2,1)/a(1,1)*qr(1,i1)))**2+& -& (f_coeff(3,1,1,1,iw)*(qr(1,i1)+a(2,1)/a(1,1)*qr(2,i1)))**2)) + else ! - !Evaluate W - RIM_acc = RIM_acc + rfac*vslab2*func*vslab2/(r1**2*(1-vslab2*func)) - !Anisotropy contribution - anis_fact=dot_product(em1_anis,(2._SP*pi*qr(:,i1)/(alat(:)*r1))**2) - func = func*anis_fact - RIM_acc_anis = RIM_acc_anis + rfac*vslab2*func*vslab2/(r1**2*(1-vslab2*func)) + ! Compute at finite frequency ! - !Store the square root of vslab - vslab(i1,1) = sqrt(vslab2)/r1 + select case (trim(rimw_type)) ! - enddo + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("semiconductor","Dirac") + ! GS: Identical to static case except for anistropy here missing, + ! merge of the two cases? + do i1=1,RIM_n_rand_pts + ! + call c2a(v_in=qr(:,i1),v_out=qr_rlu(:),mode='ki2a') + ! + slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) + !kxy + r1=sqrt((qr(1,i1)*2.*pi/alat(1))**2+(qr(2,i1)*2.*pi/alat(2))**2 + & +& (qr(3,i1)*2.*pi/alat(3))**2) + ! + !Regularization + if (r1 < 1.e-5) then + if (cut_is_slab) then + ! + ! System 2D, prefactor of slab truncation + ! + vslab(i1,1) = sqrt(4._DP*pi*(1.-exp(-q0_def_norm*lcut))/q0_def_norm**2) + RIM_acc = RIM_acc + rfac*f_coeff(1,1,1,1,iw)*(4*pi*lcut)**2 + ! + else + ! + ! System 3D, standard bare coulomb + ! + vslab(i1,1) = sqrt(4._DP*pi/q0_def_norm**2) + RIM_acc = RIM_acc + rfac*f_coeff(1,1,1,1,iw)/(1.0-(4._DP*pi)* & +& f_coeff(1,1,1,1,iw))*(4*pi)**2/q0_def_norm**2 + ! + endif + ! + cycle + ! + endif + ! + !Evaluate v_slab and interpolation function + ! + if (cut_is_slab) then + ! + vslab2 = 4._DP*pi*(1.-exp(-slab_vplane1*lcut)) + func = f_coeff(1,1,1,1,iw)*exp(sign(1._DP,real(f_coeff(idir(2)+1,1,1,1,iw),DP)) & +& *sqrt((f_coeff(idir(2)+1,1,1,1,iw)*qr_rlu(idir(2)))**2 & +& +(f_coeff(idir(3)+1,1,1,1,iw)*qr_rlu(idir(3)))**2)) + ! + else + ! + vslab2 = 4._DP*pi + func = f_coeff(1,1,1,1,iw)*exp(sign(1._DP,real(f_coeff(2,1,1,1,iw),DP)) & +& *sqrt((f_coeff(2,1,1,1,iw)*qr_rlu(1))**2 + (f_coeff(3,1,1,1,iw) & +& *qr_rlu(2))**2 + (f_coeff(4,1,1,1,iw)*qr_rlu(3))**2)) + ! + endif + ! + !Evaluate W + RIM_acc = RIM_acc + rfac*vslab2*func*vslab2/(1-vslab2*func)/r1**2 + ! + !Store the square root of vslab + vslab(i1,1) = sqrt(vslab2)/r1 + ! + enddo + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("metal") + ! + call c2a(v_in=qr(:,i1),v_out=qr_rlu(:),mode='ki2a') + ! + do i1=1,RIM_n_rand_pts + ! + r1=sqrt((qr(1,i1)*2.*pi/alat(1))**2+(qr(2,i1)*2.*pi/alat(2))**2 + & +& (qr(3,i1)*2.*pi/alat(3))**2) + ! + slab_vplane1=sqrt((qr(idir(2),i1)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i1)*2.*pi/alat(idir(3)))**2) + ! + func = f_coeff(1,1,1,1,iw) + ! + if (r1 < 1.e-5) then + ! + if (cut_is_slab) then + ! + ! System 2D, prefactor of slab truncation + ! + vslab(i1,1) = sqrt(4._DP*pi*(1.0_DP-exp(-q0_def_norm*lcut))/q0_def_norm**2) + RIM_acc = RIM_acc + rfac*func*(4.0_DP*pi*lcut)**2 + ! + else + ! + ! System 3D, standard bare coulomb + ! + vslab(i1,1) = sqrt(4._DP*pi/q0_def_norm**2) + RIM_acc = RIM_acc + rfac*func/(1.0_DP-(4._DP*pi)*func)*(4*pi)**2/q0_def_norm**2 + ! + endif + ! + cycle + ! + endif + ! + vslab2 = 4._DP*pi + ! + if (cut_is_slab) then + vslab2 = vslab2*(1.-exp(-slab_vplane1*lcut)) + func = (func+f_coeff(idir(2)+4,1,1,1,iw)*abs(qr_rlu(idir(2)))+f_coeff(idir(3)+4,1,1,1,iw) & +& *abs(qr_rlu(idir(3))))*exp(sign(1._DP,real(f_coeff(idir(2)+1,1,1,1,iw),DP)) & +& *sqrt((f_coeff(idir(2)+1,1,1,1,iw)*qr_rlu(idir(2)))**2+(f_coeff(idir(3)+1,1,1,1,iw) & +& *qr_rlu(idir(3)))**2)) + ! + else + ! + func = (func+f_coeff(5,1,1,1,iw)*abs(qr_rlu(1)) + f_coeff(6,1,1,1,iw)*abs(qr_rlu(2)) & +& + f_coeff(7,1,1,1,iw)*abs(qr_rlu(3)))*exp(sign(1._DP,real(f_coeff(3,1,1,1,iw),DP)) & +& * sqrt((f_coeff(2,1,1,1,iw)*qr_rlu(1))**2 + (f_coeff(3,1,1,1,iw)*qr_rlu(2))**2 & +& + (f_coeff(4,1,1,1,iw)*qr_rlu(3))**2)) + ! + endif + ! + !Evaluate W + RIM_acc = RIM_acc + rfac*vslab2*func*vslab2/(1-vslab2*func)/r1**2 + ! + !Store the square root of vslab + vslab(i1,1) = sqrt(vslab2)/r1 + ! + enddo + ! + end select + ! + endif ! if (RIM_id_epsm1_reference == 0) then RIM_W(iw,1,1,1)=RIM_acc else RIM_W(iw,1,1,1)=RIM_acc_anis - end if + endif + ! enddo ! + !if (RIM_W_is_diagonal) return + ! ! wings q == 0 ! - if (RIM_W_is_diagonal) return do iw=1,Xw%n_freqs do i2=2,RIM_W_ng + ! + RIM_acc=0._DP + ! + do i3=1,RIM_n_rand_pts + ! + call c2a(v_in=qr(:,i3),v_out=qr_rlu(:),mode='ki2a') + ! + r1=sqrt((qr(1,i3)*2.*pi/alat(1))**2+(qr(2,i3)*2.*pi/alat(2))**2 + & +& (qr(3,i3)*2.*pi/alat(3))**2) + ! + if ( r1 < 1e-5) then + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! Small q behaviour differs for metals and semiconductors + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + select case (trim(rimw_type)) + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("metal") + ! + if (cut_is_slab) then + slab_vz1=g_vec(i2,idir(1))*2.*pi/alat(idir(1)) + slab_vplane1=sqrt((g_vec(i2,idir(2))*2.*pi/alat(idir(2)))**2 + & +& (g_vec(i2,idir(3))*2.*pi/alat(idir(3)))**2) + if (slab_vplane1 < q0_def_norm) slab_vplane1 = q0_def_norm + vslab2 = 4._DP*pi*(1.-exp(-slab_vplane1*lcut)*cos(slab_vz1*lcut))**0.5 + ! + do j=1,2 + ! Load q_grid_b in iku + qtmp(:)=q_grid_b_iku(idir(1+j),:) + ! + if (iw ==1) then + ! System 2D, real part of f func behaves as fcoeff*q, if imag present as vec(q) + func = (real(f_coeff(2,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),1)) + & +& real(f_coeff(3,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),2)) + & +& real(f_coeff(4,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),3)) + & +& aimag(f_coeff(2,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),1) + & +& aimag(f_coeff(3,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),2) + & +& aimag(f_coeff(4,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),3)) & +& /iku_v_norm(qtmp(:)) + ! + RIM_acc = RIM_acc + rfac*vslab2**2*func/(slab_vplane1**2+slab_vz1**2)/2.0_SP/q0_def_norm + ! + else + if ((g_vec(i2,idir(2))==g_vec(i2,idir(3))) .and. (g_vec(i2,idir(2))==0)) then + func = ((f_coeff(2,1,i2,iq,iw)-f_coeff(5,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),1)**2 + & +& (f_coeff(3,1,i2,iq,iw)-f_coeff(6,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),2)**2 + & +& (f_coeff(4,1,i2,iq,iw)-f_coeff(7,1,i2,iq,iw))*q_grid_b_rlu(idir(1+j),3)**2 ) & +& /iku_v_norm(qtmp(:))**2/2.0_SP + ! + ! AF: is this correct ?? + ! GS: yes, the out of plane wings do not contribute. In hindsight, this part is superfluous + ! + RIM_acc = RIM_acc + ! + else + func = ((f_coeff(2,1,i2,iq,iw)-f_coeff(5,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),1)) + & +& (f_coeff(3,1,i2,iq,iw)-f_coeff(6,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),2)) + & +& (f_coeff(4,1,i2,iq,iw)-f_coeff(7,1,i2,iq,iw))*abs(q_grid_b_rlu(idir(1+j),3)) ) & +& /iku_v_norm(qtmp(:))/2.0_SP + ! + RIM_acc = RIM_acc + rfac*vslab2**2*func/(slab_vplane1**2+slab_vz1**2) & +& /2.0_SP/q0_def_norm + endif + endif + ! + enddo + ! + else + ! + ! System 3D, f func behaves as f = fcoeff*q^2 + ! we loop on 3D directions as fcoeff may be anisotropic + ! + do j=1,3 + ! Load q_grid_b in iku + qtmp(:)=q_grid_b_iku(j,:) + ! + func = (qtmp(1)*(qtmp(1)*f_coeff(5,1,i2,iq,iw) + 2_SP*qtmp(2)* & +& f_coeff(6,1,i2,iq,iw)) + qtmp(2)*(qtmp(2)*f_coeff(8,1,i2,iq,iw) & +& + 2_SP*qtmp(3)*f_coeff(9,1,i2,iq,iw)) +qtmp(3)*(qtmp(3) & +& *f_coeff(10,1,i2,iq,iw) + 2_SP*qtmp(1)*f_coeff(7,1,i2,iq,iw))) & +& /iku_v_norm(qtmp(:))**2 + ! + slab_vplane1=sqrt((g_vec(i2,1)*2.0_SP*pi/alat(1))**2 + (g_vec(i2,2)*2.0_SP*pi & +& /alat(2))**2 + (g_vec(i2,3)*2.*pi/alat(3))**2) + ! + RIM_acc = RIM_acc + rfac*(4._SP*pi)**2*func/slab_vplane1**2/3.0 + ! + enddo + ! + endif + ! + cycle + ! + end select + ! + endif + ! + !Evaluate func + ! + func = f_coeff(1,1,i2,iq,iw) & +& + qr(1,i3)*(f_coeff(2,1,i2,iq,iw)+qr(1,i3)*f_coeff(5,1,i2,iq,iw) & +& + 2._DP*qr(2,i3)*f_coeff(6,1,i2,iq,iw)) & +& + qr(2,i3)*(f_coeff(3,1,i2,iq,iw)+qr(2,i3)*f_coeff(8,1,i2,iq,iw) & +& + 2._DP*qr(3,i3)*f_coeff(9,1,i2,iq,iw)) & +& + qr(3,i3)*(f_coeff(4,1,i2,iq,iw)+qr(3,i3)*f_coeff(10,1,i2,iq,iw) & +& + 2._DP*qr(1,i3)*f_coeff(7,1,i2,iq,iw)) & +& +qr(1,i3)**2*(qr(1,i3)*f_coeff(11,1,i2,iq,iw)+3_DP*qr(2,i3)* & +& f_coeff(12,1,i2,iq,iw)+3_DP*qr(3,i3)*f_coeff(13,1,i2,iq,iw)) & +& +qr(2,i3)**2*(qr(2,i3)*f_coeff(15,1,i2,iq,iw)+3_DP*qr(3,i3)* & +& f_coeff(16,1,i2,iq,iw)+3_DP*qr(1,i3)*f_coeff(14,1,i2,iq,iw)) & +& +qr(3,i3)**2*(qr(3,i3)*f_coeff(19,1,i2,iq,iw)+3_DP*qr(1,i3)* & +& f_coeff(17,1,i2,iq,iw)+3_DP*qr(2,i3)*f_coeff(18,1,i2,iq,iw)) & +& +6*qr(1,i3)*qr(2,i3)*qr(3,i3)*f_coeff(20,1,i2,iq,iw) + ! + select case (trim(rimw_type)) + ! + case("metal") + ! + if (cut_is_slab) then + ! + if (iw ==1) then + func = real(f_coeff(2,1,i2,iq,iw))*abs(qr_rlu(1)) + real(f_coeff(3,1,i2,iq,iw))*abs(qr_rlu(2)) & +& + real(f_coeff(4,1,i2,iq,iw))*abs(qr_rlu(3)) + aimag(f_coeff(2,1,i2,iq,iw))*qr_rlu(1) & +& + aimag(f_coeff(3,1,i2,iq,iw))*qr_rlu(2) + aimag(f_coeff(4,1,i2,iq,iw))*qr_rlu(3) & +& + f_coeff(5,i2,1,iq,iw)*(qr_rlu(1))**2 + f_coeff(6,i2,1,iq,iw)*(qr_rlu(2))**2 & +& + f_coeff(7,i2,1,iq,iw)*(qr_rlu(3))**2 + f_coeff(8,i2,1,iq,iw) & +& * sign(qr_rlu(1)**2,qr_rlu(1)) + f_coeff(9,i2,1,iq,iw)*sign(qr_rlu(2)**2,qr_rlu(2)) & +& + f_coeff(10,i2,1,iq,iw)*sign(qr_rlu(3)**2,qr_rlu(3)) + else + ! + if ( (g_vec(i2,idir(2))==g_vec(i2,idir(3))) .and. (g_vec(i2,idir(2))==0)) then + func = ((f_coeff(2,1,i2,iq,iw)+f_coeff(5,1,i2,iq,iw)+ (f_coeff(2,1,i2,iq,iw) & +& -f_coeff(5,1,i2,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**2 + & +& (f_coeff(3,1,i2,iq,iw)+f_coeff(6,1,i2,iq,iw)+ (f_coeff(3,1,i2,iq,iw) & +& -f_coeff(6,1,i2,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**2 + & +& (f_coeff(4,1,i2,iq,iw)+f_coeff(7,1,i2,iq,iw)+ (f_coeff(4,1,i2,iq,iw) & +& -f_coeff(7,1,i2,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**2 + & +& (f_coeff(8,1,i2,iq,iw)+f_coeff(11,1,i2,iq,iw)+ (f_coeff(8,1,i2,iq,iw) & +& -f_coeff(11,1,i2,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**3 + & +& (f_coeff(9,1,i2,iq,iw)+f_coeff(12,1,i2,iq,iw)+ (f_coeff(9,1,i2,iq,iw) & +& -f_coeff(12,1,i2,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**3 + & +& (f_coeff(10,1,i2,iq,iw)+f_coeff(13,1,i2,iq,iw)+ (f_coeff(10,1,i2,iq,iw) & +& -f_coeff(13,1,i2,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**3)/2.0_SP + else + ! + func = ((f_coeff(2,1,i2,iq,iw)+f_coeff(5,1,i2,iq,iw)+ (f_coeff(2,1,i2,iq,iw) & +& -f_coeff(5,1,i2,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1) + & +& (f_coeff(3,1,i2,iq,iw)+f_coeff(6,1,i2,iq,iw)+ (f_coeff(3,1,i2,iq,iw) & +& -f_coeff(6,1,i2,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2) + & +& (f_coeff(4,1,i2,iq,iw)+f_coeff(7,1,i2,iq,iw)+ (f_coeff(4,1,i2,iq,iw) & +& -f_coeff(7,1,i2,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3) + & +& (f_coeff(8,1,i2,iq,iw)+f_coeff(11,1,i2,iq,iw)+ (f_coeff(8,1,i2,iq,iw) & +& -f_coeff(11,1,i2,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**2 + & +& (f_coeff(9,1,i2,iq,iw)+f_coeff(12,1,i2,iq,iw)+ (f_coeff(9,1,i2,iq,iw) & +& -f_coeff(12,1,i2,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**2 + & +& (f_coeff(10,1,i2,iq,iw)+f_coeff(13,1,i2,iq,iw)+ (f_coeff(10,1,i2,iq,iw) & +& -f_coeff(13,1,i2,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**2)/2.0_SP + endif + ! + endif + ! + endif + ! + end select + ! + r1=sqrt((qr(1,i3)*2.*pi/alat(1))**2+(qr(2,i3)*2.*pi/alat(2))**2 + & +& (qr(3,i3)*2.*pi/alat(3))**2) + slab_vplane1=sqrt((qr(idir(2),i3)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i3)*2.*pi/alat(idir(3)))**2) + ! + if (r1 < 1.e-5) then + ! + vslab(i3,1) = sqrt(4._DP*pi/q0_def_norm**2) + ! + if (cut_is_slab) vslab(i3,1) = vslab(i3,1)*sqrt(1.-exp(-q0_def_norm*lcut)) + ! + else + ! + vslab(i3,1) =sqrt(4._DP*pi/r1**2) + ! + if (cut_is_slab) vslab(i3,1) = vslab(i3,1)*sqrt(1.-exp(-slab_vplane1*lcut)) + ! + endif + ! + ! Accumulate Wc + ! + RIM_acc = RIM_acc + rfac*vslab(i3,1)**2*func*vslab(i3,i2)**2/(1._DP- & +& vslab(i3,1)*func*vslab(i3,i2)) + ! + enddo + RIM_W(iw,iq,1,i2)=RIM_acc + ! old version symmetrization RIM_W(iw,iq,i2,1)=RIM_W(iw,iq,1,i2) + enddo + ! + enddo + ! + ! + do iw=1,Xw%n_freqs + do i2=2,RIM_W_ng + ! RIM_acc=0._DP ! do i3=1,RIM_n_rand_pts ! - !Evaluate func + call c2a(v_in=qr(:,i3),v_out=qr_rlu(:),mode='ki2a') + ! + r1=sqrt((qr(1,i3)*2.*pi/alat(1))**2+(qr(2,i3)*2.*pi/alat(2))**2 + & +& (qr(3,i3)*2.*pi/alat(3))**2) ! - func = f_coeff(1,1,i2,iq,iw)+qr(1,i3)*(f_coeff(2,1,i2,iq,iw)+qr(1,i3)*f_coeff(4,1,i2,iq,iw))& -& +qr(2,i3)*(f_coeff(3,1,i2,iq,iw)+qr(2,i3)*f_coeff(6,1,i2,iq,iw)& -& +2._DP*qr(1,i3)*f_coeff(5,1,i2,iq,iw)) ! - RIM_acc = RIM_acc + rfac*vslab(i3,1)**2*func*vslab(i3,i2)**2/(1._DP-vslab(i3,1)*func*vslab(i3,i2)) + if ( r1 < 1e-5) then + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! Small q behaviour differs for metals and semiconductors + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + select case (trim(rimw_type)) + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + case("metal") + ! + if (cut_is_slab) then + ! + slab_vz1=g_vec(i2,idir(1))*2.*pi/alat(idir(1)) + slab_vplane1=sqrt((g_vec(i2,idir(2))*2.*pi/alat(idir(2)))**2 + & +& (g_vec(i2,idir(3))*2.*pi/alat(idir(3)))**2) + if (slab_vplane1 < q0_def_norm) slab_vplane1 = q0_def_norm + vslab2 = 4._DP*pi*(1.-exp(-slab_vplane1*lcut)*cos(slab_vz1*lcut))**0.5 + ! + do j=1,2 + ! Load q_grid_b in iku + qtmp(:)=q_grid_b_iku(idir(1+j),:) + ! + if (iw ==1) then + ! System 2D, real part of f func behaves as fcoeff*q, if imag present as vec(q) + func = (real(f_coeff(2,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),1)) + & +& real(f_coeff(3,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),2)) + & +& real(f_coeff(4,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),3)) + & +& aimag(f_coeff(2,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),1) + & +& aimag(f_coeff(3,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),2) + & +& aimag(f_coeff(4,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),3)) & +& /iku_v_norm(qtmp(:)) + ! + RIM_acc = RIM_acc + rfac*vslab2**2*func/(slab_vplane1**2+slab_vz1**2)/2.0_SP/q0_def_norm + ! + else + ! + if ((g_vec(i2,idir(2))==g_vec(i2,idir(3))) .and. (g_vec(i2,idir(2))==0)) then + func = ( (f_coeff(2,i2,1,iq,iw)+f_coeff(5,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),1)**2 + & +& (f_coeff(3,i2,1,iq,iw)+f_coeff(6,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),2)**2 + & +& (f_coeff(4,i2,1,iq,iw)+f_coeff(7,i2,1,iq,iw))*q_grid_b_rlu(idir(1+j),3)**2 ) & +& /iku_v_norm(qtmp(:))**2/2.0_SP + ! + ! AF: is this correct ? + ! GS: yes, the out of plane wings do not contribute. In hindsight, this part is superfluous + ! + RIM_acc = RIM_acc + ! + else + ! + func = ( (f_coeff(2,i2,1,iq,iw)-f_coeff(5,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),1)) + & +& (f_coeff(3,i2,1,iq,iw)-f_coeff(6,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),2)) + & +& (f_coeff(4,i2,1,iq,iw)-f_coeff(7,i2,1,iq,iw))*abs(q_grid_b_rlu(idir(1+j),3)) ) & +& /iku_v_norm(qtmp(:))/2.0_SP + ! + RIM_acc = RIM_acc + rfac*vslab2**2*func/(slab_vplane1**2+slab_vz1**2) & +& /2.0_SP/q0_def_norm + endif + ! + endif + ! + enddo + ! + else + ! + ! System 3D, f func behaves as f = fcoeff*q^2 + ! we loop on 3D directions as fcoeff may be anisotropic + ! + do j=1,3 + ! Load q_grid_b in iku + qtmp(:)=q_grid_b_iku(j,:) + ! + func = (qtmp(1)*(qtmp(1)*f_coeff(5,i2,1,iq,iw) + 2_SP*qtmp(2)* & +& f_coeff(6,i2,1,iq,iw)) + qtmp(2)*(qtmp(2)*f_coeff(8,i2,1,iq,iw) & +& + 2_SP*qtmp(3)*f_coeff(9,i2,1,iq,iw)) +qtmp(3)*(qtmp(3) & +& *f_coeff(10,i2,1,iq,iw) + 2_SP*qtmp(1)*f_coeff(7,i2,1,iq,iw))) & +& /iku_v_norm(qtmp(:))**2 + ! + slab_vplane1=sqrt((g_vec(i2,1)*2.0_SP*pi/alat(1))**2 + (g_vec(i2,2)*2.0_SP*pi & +& /alat(2))**2 + (g_vec(i2,3)*2.*pi/alat(3))**2) + ! + RIM_acc = RIM_acc + rfac*(4._SP*pi)**2*func/slab_vplane1**2/3.0 + ! + enddo + ! + endif + cycle + ! + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + !case("semiconductor") (removed specific case for now, computed with general rule) + ! + !Load q_grid_b in rlu + !if (cut_is_slab) then + ! + ! System 2D, f func behaves as f = fcoeff*sqrt(q) + ! fcoeff may differ in the planar directions + ! + !do j=1,2 + ! + !qr(:,i3)=q_grid_b_iku(idir(1+j),:) + ! + !func = (qr(1,i3)*f_coeff(2,1,i2,iq,iw) + qr(2,i3)*f_coeff(3,1,i2,iq,iw) + ! + qr(3,i3)*f_coeff(4,1,i2,iq,iw))/sqrt(iku_v_norm(qr(:,i3))) + ! + !slab_vz1=g_vec(i2,idir(1))*2.*pi/alat(idir(1)) + !slab_vplane1=sqrt(g_vec(i2,idir(2))*2.*pi/alat(idir(2))**2 + + !g_vec(i2,idir(3))*2.*pi/alat(idir(3))**2) + !RIM_acc = RIM_acc + rfac*(4._SP*pi)**2*func/(slab_vplane1**2+slab_vz1**2)* + !(1-exp(-lcut*slab_vplane1)*cos(lcut*slab_vz1))**2/2.0_SP + ! + !end do + !else + ! + ! System 3D, f func behaves as f = fcoeff*q + ! fcoeff may differ in the three directions + ! + !do j=1,3 + ! + !qr(:,i3)=q_grid_b_iku(j,:) + ! + !func = (qr(1,i3)*f_coeff(2,1,i2,iq,iw) + qr(2,i3)*f_coeff(3,1,i2,iq,iw) + ! + qr(3,i3)*f_coeff(4,1,i2,iq,iw))/iku_v_norm(qr(:,i3)) + ! + !slab_vplane1=sqrt(g_vec(i2,idir(1))*2.*pi/alat(idir(1))**2 + + !g_vec(i2,idir(2))*2.*pi/alat(idir(2))**2 +g_vec(i2,idir(3))*2. + ! *pi/alat(idir(3))**2) + !RIM_acc = RIM_acc + rfac*(4._SP*pi)**2*func/slab_vplane1/q0_def_norm/3.0 + ! + ! + !end do + !end if + ! + end select + ! + endif + ! + !Evaluate func + ! + func = f_coeff(1,i2,1,iq,iw) & +& + qr(1,i3)*(f_coeff(2,i2,1,iq,iw)+qr(1,i3)*f_coeff(5,i2,1,iq,iw) & +& + 2._DP*qr(2,i3)*f_coeff(6,i2,1,iq,iw)) & +& + qr(2,i3)*(f_coeff(3,i2,1,iq,iw)+qr(2,i3)*f_coeff(8,i2,1,iq,iw) & +& + 2._DP*qr(3,i3)*f_coeff(9,i2,1,iq,iw)) & +& + qr(3,i3)*(f_coeff(4,i2,1,iq,iw)+qr(3,i3)*f_coeff(10,i2,1,iq,iw) & +& + 2._DP*qr(1,i3)*f_coeff(7,i2,1,iq,iw)) & +& +qr(1,i3)**2*(qr(1,i3)*f_coeff(11,i2,1,iq,iw)+3_DP*qr(2,i3)* & +& f_coeff(12,i2,1,iq,iw)+3_DP*qr(3,i3)*f_coeff(13,i2,1,iq,iw)) & +& +qr(2,i3)**2*(qr(2,i3)*f_coeff(15,i2,1,iq,iw)+3_DP*qr(3,i3)* & +& f_coeff(16,i2,1,iq,iw)+3_DP*qr(1,i3)*f_coeff(14,i2,1,iq,iw)) & +& +qr(3,i3)**2*(qr(3,i3)*f_coeff(19,i2,1,iq,iw)+3_DP*qr(1,i3)* & +& f_coeff(17,i2,1,iq,iw)+3_DP*qr(2,i3)*f_coeff(18,i2,1,iq,iw)) & +& +6*qr(1,i3)*qr(2,i3)*qr(3,i3)*f_coeff(20,i2,1,iq,iw) + ! + select case (trim(rimw_type)) + ! + case("metal") + ! + if (cut_is_slab) then + ! + if (iw ==1) then + func = real(f_coeff(2,i2,1,iq,iw))*abs(qr_rlu(1)) + real(f_coeff(3,i2,1,iq,iw))*abs(qr_rlu(2)) & +& + real(f_coeff(4,i2,1,iq,iw))*abs(qr_rlu(3)) + aimag(f_coeff(2,i2,1,iq,iw))*qr_rlu(1) & +& + aimag(f_coeff(3,i2,1,iq,iw))*qr_rlu(2) + aimag(f_coeff(4,i2,1,iq,iw))*qr_rlu(3) & +& + f_coeff(5,i2,1,iq,iw)*(qr_rlu(1))**2 + f_coeff(6,i2,1,iq,iw)*(qr_rlu(2))**2 & +& + f_coeff(7,i2,1,iq,iw)*(qr_rlu(3))**2 + f_coeff(8,i2,1,iq,iw) & +& * sign(qr_rlu(1)**2,qr_rlu(1)) + f_coeff(9,i2,1,iq,iw)*sign(qr_rlu(2)**2,qr_rlu(2)) & +& + f_coeff(10,i2,1,iq,iw)*sign(qr_rlu(3)**2,qr_rlu(3)) + else + ! + if ( (g_vec(i2,idir(2))==g_vec(i2,idir(3))) .and. (g_vec(i2,idir(2))==0)) then + func = ((f_coeff(2,i2,1,iq,iw)+f_coeff(5,i2,1,iq,iw)+ (f_coeff(2,i2,1,iq,iw) & +& -f_coeff(5,i2,1,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**2 + & +& (f_coeff(3,i2,1,iq,iw)+f_coeff(6,i2,1,iq,iw)+ (f_coeff(3,i2,1,iq,iw) & +& -f_coeff(6,i2,1,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**2 + & +& (f_coeff(4,i2,1,iq,iw)+f_coeff(7,i2,1,iq,iw)+ (f_coeff(4,i2,1,iq,iw) & +& -f_coeff(7,i2,1,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**2 + & +& (f_coeff(8,i2,1,iq,iw)+f_coeff(11,i2,1,iq,iw)+ (f_coeff(8,i2,1,iq,iw) & +& -f_coeff(11,i2,1,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**3 + & +& (f_coeff(9,i2,1,iq,iw)+f_coeff(12,i2,1,iq,iw)+ (f_coeff(9,i2,1,iq,iw) & +& -f_coeff(12,i2,1,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**3 + & +& (f_coeff(10,i2,1,iq,iw)+f_coeff(13,i2,1,iq,iw)+ (f_coeff(10,i2,1,iq,iw) & +& -f_coeff(13,i2,1,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**3)/2.0_SP + else + ! + func = ((f_coeff(2,i2,1,iq,iw)+f_coeff(5,i2,1,iq,iw)+ (f_coeff(2,i2,1,iq,iw) & +& -f_coeff(5,i2,1,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1) + & +& (f_coeff(3,i2,1,iq,iw)+f_coeff(6,i2,1,iq,iw)+ (f_coeff(3,i2,1,iq,iw) & +& -f_coeff(6,i2,1,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2) + & +& (f_coeff(4,i2,1,iq,iw)+f_coeff(7,i2,1,iq,iw)+ (f_coeff(4,i2,1,iq,iw) & +& -f_coeff(7,i2,1,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3) + & +& (f_coeff(8,i2,1,iq,iw)+f_coeff(11,i2,1,iq,iw)+ (f_coeff(8,i2,1,iq,iw) & +& -f_coeff(11,i2,1,iq,iw))*sign(1._SP,qr_rlu(1)))*qr_rlu(1)**2 + & +& (f_coeff(9,i2,1,iq,iw)+f_coeff(12,i2,1,iq,iw)+ (f_coeff(9,i2,1,iq,iw) & +& -f_coeff(12,i2,1,iq,iw))*sign(1._SP,qr_rlu(2)))*qr_rlu(2)**2 + & +& (f_coeff(10,i2,1,iq,iw)+f_coeff(13,i2,1,iq,iw)+ (f_coeff(10,i2,1,iq,iw) & +& -f_coeff(13,i2,1,iq,iw))*sign(1._SP,qr_rlu(3)))*qr_rlu(3)**2)/2.0_SP + endif + ! + endif + ! + endif + ! + end select + ! + r1=sqrt((qr(1,i3)*2.*pi/alat(1))**2+(qr(2,i3)*2.*pi/alat(2))**2 + & +& (qr(3,i3)*2.*pi/alat(3))**2) + slab_vplane1=sqrt((qr(idir(2),i3)*2.*pi/alat(idir(2)))**2 + & +& (qr(idir(3),i3)*2.*pi/alat(idir(3)))**2) + ! + if (r1 < 1.e-5) then + ! + vslab(i3,1) = sqrt(4._DP*pi/q0_def_norm**2) + ! + if (cut_is_slab) vslab(i3,1) = vslab(i3,1)*sqrt(1.-exp(-q0_def_norm*lcut)) + ! + else + ! + vslab(i3,1) =sqrt(4._DP*pi/r1**2) + ! + if (cut_is_slab) vslab(i3,1) = vslab(i3,1)*sqrt(1.-exp(-slab_vplane1*lcut)) + ! + endif + ! + ! Accumulate Wc + ! + RIM_acc = RIM_acc + rfac*vslab(i3,1)**2*func*vslab(i3,i2)**2/(1._DP- & +& vslab(i3,1)*func*vslab(i3,i2)) + ! enddo - RIM_W(iw,iq,1,i2)=RIM_acc - RIM_W(iw,iq,i2,1)=RIM_W(iw,iq,1,i2) + ! + RIM_W(iw,iq,i2,1)=RIM_acc + ! old version symmetrization RIM_W(iw,iq,i2,1)=RIM_W(iw,iq,1,i2) enddo enddo ! + YAMBO_FREE(vslab) + ! end subroutine diff --git a/src/interface/INIT_activate.F b/src/interface/INIT_activate.F index a9c046beb0..82a3544379 100644 --- a/src/interface/INIT_activate.F +++ b/src/interface/INIT_activate.F @@ -145,9 +145,7 @@ subroutine INIT_activate() ! if (l_rim) call initactivate(1,'RandQpts RandGvec QpgFull Em1Anys IDEm1Ref') ! - !RIM_W - ! - if (l_rim_w) call initactivate(1,'RandQpts RandGvecW rimw_type') + if (l_rim_w) call initactivate(1,'RandQpts RandGvecW rimw_type Dirac_v_f rimw_check') ! !Col CUTOFF ! @@ -323,7 +321,7 @@ subroutine INIT_activate() ! if ( (l_em1d.and.l_ppa) .or. (l_em1d.and.l_mpa) .or. (l_em1s.and.l_cohsex)) then call INIT_QP_ctl_switch('X') - call initactivate(1,'RIM_W RIM_W_diag RIM_W_graph') + call initactivate(1,'RIM_W RIM_W_diag') ! #if defined _RT call INIT_RT_ctl_switch('X') @@ -370,9 +368,9 @@ subroutine INIT_activate() #if defined _QED if (l_elphoton_corr) then if (l_gw0) then - call initactivate(1,'FFTGvecs RandQpts QEDRLvcs GbndRnge GDamping dScStep DysSolver') + call initactivate(1,'FFTGvecs RandQpts QEDRLvcs GbndRnge GDamping G_mprop dScStep DysSolver') if (trim(QP_solver)=="g") then - call initactivate(1,'GEnSteps GEnRnge GEnMode GDmRnge GreenFTresh GreenF2QP') + call initactivate(1,'GEnSteps GEnRnge GEnMode GDmRnge GreenFTresh GreenF2QP') else call initactivate(1,'SCEtresh') if (.not.l_cohsex) call initactivate(1,'NewtDchk ExtendOut OnMassShell QPExpand') @@ -397,7 +395,7 @@ subroutine INIT_activate() ! if (l_gw0) then if (.not.l_cohsex.or.COHSEX_use_empties) call initactivate(1,'GbndRnge') - if (.not.l_cohsex.and.trim(QP_solver)/='g') call initactivate(1,'GDamping') + if (.not.l_cohsex.and.trim(QP_solver)/='g') call initactivate(1,'GDamping G_mprop') if (.not.l_cohsex) call initactivate(1,'dScStep') if (.not.l_elphoton_corr) then if (.not.l_ppa.and..not.l_mpa.and..not.l_cohsex) & @@ -424,16 +422,6 @@ subroutine INIT_activate() endif endif ! -#if defined _PHEL - ! - if (l_phel_corr) then - call initactivate(1,'ElecTemp BoseTemp PH_SE_mode GphBRnge GDamping PHELQpts ElPhModes PHDbGdsize DbGdQsize GDamping') - call initactivate(1,'ExtendOut OnMassShell') - if ( l_gw0) call initactivate(1,'DysSolver GEnSteps GEnRnge PHEL_QPH_En') - if (.not. l_gw0) call initactivate(1,'GDamping PHELTrans') - endif - ! -#endif #if defined _ELPH ! if (l_elph_corr) then diff --git a/src/interface/INIT_barriers.F b/src/interface/INIT_barriers.F index 703e444bac..779a8f7660 100644 --- a/src/interface/INIT_barriers.F +++ b/src/interface/INIT_barriers.F @@ -104,7 +104,7 @@ subroutine INIT_barriers() endif ! if (gw0_raxis) then - on_runlevels='dyson gw0 em1d el_el_corr el_ph_corr ph_el_corr el_photon_corr HF_and_locXC' + on_runlevels='dyson gw0 em1d el_el_corr el_ph_corr ph_el_corr el_photon_corr HF_and_locXC rim_w' call switch_off_runlevel('all',except=trim(on_runlevels)//' '//trim(always_on_runlevels)) goto 1 endif @@ -115,7 +115,7 @@ subroutine INIT_barriers() goto 1 endif if (gw0_mpa) then - on_runlevels='dyson gw0 mpa em1d el_el_corr el_ph_corr HF_and_locXC' + on_runlevels='dyson gw0 mpa em1d el_el_corr el_ph_corr HF_and_locXC rim_w' call switch_off_runlevel('all',except=trim(on_runlevels)//' '//trim(always_on_runlevels)) goto 1 endif diff --git a/src/interface/INIT_load.F b/src/interface/INIT_load.F index c2fb198625..92f3eb1a17 100644 --- a/src/interface/INIT_load.F +++ b/src/interface/INIT_load.F @@ -15,21 +15,21 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) use it_tools, ONLY:it use it_m, ONLY:initdefs,E_unit,G_unit,T_unit,Bfield_unit,MEM_unit,& & Time_unit,I_unit,Angle_unit,V_parallel,initmode,V_ph,& -& V_RL,V_kpt,V_sc,V_qp,V_io,V_general,V_resp,V_real_time,V_nl_optics +& V_RL,V_kpt,V_sc,V_qp,V_coul,V_io,V_general,V_resp,V_real_time,V_nl_optics use X_m, ONLY:Chi_mode,X_t,q_plus_G_direction,Q_Shift_Order,& & global_gauge,Chi_linalg_mode, & & X_terminator_Kind,X_terminator_E,X_DbGd_percentual use DIPOLES, ONLY:DIPOLE_t use com, ONLY:grid_paths use stderr, ONLY:slash - use QP_m, ONLY:QP_cg_percent,QP_G_damp,QP_solver,& + use QP_m, ONLY:QP_cg_percent,QP_G_damp,QP_G_mprop,QP_solver,& & QP_n_G_bands,QP_ng_Sx,QP_ng_Sc,QP_ng_SH,QP_ng_Vxc,GW_terminator_E,GW_terminator_Kind,& & QP_G_er,QP_G_ir,QP_G_dr,QP_Sc_steps,QP_G_solver,& & QP_dSc_delta,QP_G_Zoom_treshold,GF_energy_range_mode use LIVE_t, ONLY:nhash use wave_func, ONLY:wf_ng use D_lattice, ONLY:Tel,non_periodic_directions,molecule_position,Bose_Temp - use R_lattice, ONLY:ng_closed,QP_states_k,nXkibz,k_GRIDS_string,RIM_W_ng,rimw_type,& + use R_lattice, ONLY:ng_closed,QP_states_k,nXkibz,k_GRIDS_string,RIM_W_ng,rimw_type,Dirac_v_f,rimw_check, & & bz_samp,RIM_ng,RIM_epsm1,RIM_id_epsm1_reference,& & RIM_n_rand_pts,cyl_ph_radius,box_length,cyl_length,cut_geometry,ws_cutoff use BS, ONLY:BSE_mode,BSE_prop,BSK_mode,BS_eh_en,BS_eh_win,BS_q,BS_bands,& @@ -258,6 +258,7 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) call it(defs,'GDamping', '[GW] G[W] damping',QP_G_damp,E_unit) #else call it(defs,'GDamping', '[GW] G[W] damping',QP_G_damp,E_unit,verb_level=V_qp) + call it(defs,'G_mprop', '[GW] G[W] propagator order',QP_G_mprop) #endif call it(defs,'GDmRnge', '[GW] G_gw damping range',QP_G_dr,E_unit) call it(defs,'dScStep', '[GW] Energy step to evaluate Z factors',QP_dSc_delta,E_unit,verb_level=V_qp) @@ -396,8 +397,12 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) call it('f',defs,'QpgFull', '[F RIM] Coulomb interaction: Full matrix',verb_level=V_RL) ! ! RIM-W - call it(defs,'RandGvecW','[RIM-W] Screened interaction RS components',RIM_W_ng,G_unit) - call it(defs,'rimw_type','[RIM-W] Screened interaction limit metal/semiconductor/graphene',rimw_type) + ! + call it(defs,'RandGvecW','[RIM-W] Screened interaction RS components (smaller than X matrix)',RIM_W_ng,G_unit) + call it(defs,'rimw_type','[RIM-W] Screened interaction limit metal/semiconductor/2D Dirac',rimw_type) + call it(defs,'Dirac_v_f','[RIM-W] In case of 2D Dirac, specify the Fermi velocity of the cone (10^6 m/s)',& +& Dirac_v_f,verb_level=V_coul) + call it('f',defs,'rimw_check','[RIM-W] Debug option to check W, XV and auxiliary functions',verb_level=V_coul) ! ! CUTOFF ! diff --git a/src/interface/INIT_read_command_line.F b/src/interface/INIT_read_command_line.F index 30ab9e718f..b1ab8b358b 100644 --- a/src/interface/INIT_read_command_line.F +++ b/src/interface/INIT_read_command_line.F @@ -19,8 +19,8 @@ subroutine INIT_read_command_line(rstr,init_) use LIVE_t, ONLY:USER_wall_time_string,GET_user_WALL_time use stderr, ONLY:STRING_split,STRING_match,STRING_same use it_tools, ONLY:runlevel_is_on,switch_off_runlevel - use it_m, ONLY:V_RL,V_kpt,V_sc,V_qp,V_io,V_general,V_resp, & -& V_real_time,V_nl_optics,V_all,V_parallel,V_ph, & + use it_m, ONLY:V_RL,V_kpt,V_sc,V_qp,V_coul,V_io,V_general,V_resp, & +& V_real_time,V_nl_optics,V_all,V_parallel,V_ph, & & infile_verbosity,nrnlvls,rnlvls #if defined _SC || defined _RT use hamiltonian, ONLY:H_potential @@ -122,6 +122,7 @@ subroutine INIT_read_command_line(rstr,init_) ! V_parallel=9 ! V_nl_optics=10 ! V_ph=11 + ! V_coul=12 ! V_all=99 ! if ( trim(rstr_piece(i1)) == 'infver' ) then @@ -136,6 +137,7 @@ subroutine INIT_read_command_line(rstr,init_) if (STRING_same( trim(rstr_piece(i1+1)) , "nl")) infile_verbosity=V_nl_optics if (STRING_same( trim(rstr_piece(i1+1)) , "ph")) infile_verbosity=V_ph if (STRING_same( trim(rstr_piece(i1+1)) , "par")) infile_verbosity=V_parallel + if (STRING_same( trim(rstr_piece(i1+1)) , "coul")) infile_verbosity=V_coul if (STRING_same( trim(rstr_piece(i1+1)) , "all")) infile_verbosity=V_all endif ! diff --git a/src/io/io_MPA.F b/src/io/io_MPA.F index f6f226284d..76ae18c96d 100644 --- a/src/io/io_MPA.F +++ b/src/io/io_MPA.F @@ -105,7 +105,7 @@ integer function io_MPA(X,Xw,ID) call def_variable_elemental(ID_frag,trim(ch),sec_size,SP,1,par_io_kind='independent') ! call io_variable_elemental(ID_frag, VAR=" :: Current Q-pt index :",I0=iq) - call io_variable_elemental(ID_frag, VAR=" :: Number of poles :",I0=Xw%n_freqs/2,CHECK=.true.,OP=(/"=="/)) + call io_variable_elemental(ID_frag, VAR=" :: Number of poles :",I0=X%mpa_npoles,CHECK=.true.,OP=(/"=="/)) ! call def_variable_elemental(ID_frag,"",0,0,1) io_MPA=io_status(ID_frag) @@ -129,35 +129,15 @@ integer function io_MPA(X,Xw,ID) ! ch="MPA_E_Q_"//trim(intc(iq)) ! - if(write_is_on(ID)) then -#if defined _PAR_IO - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP,par_io_kind='collective') -#else - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP,par_io_kind='independent') -#endif - call io_variable_bulk(ID_frag,1,C3=MPA_E_par(1)%blc(:,:,:Xw%n_freqs),IPOS=(/1,MPA_E_par(1)%rows(1),MPA_E_par(1)%cols(1),1/)) - ! - else if(read_is_on(ID)) then - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP) - call io_variable_bulk(ID_frag,1,C3=MPA_E_par(1)%blc(:,:,:Xw%n_freqs),IPOS=(/1,MPA_E_par(1)%rows(1),MPA_E_par(1)%cols(1),1/)) - endif + call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,X%mpa_npoles/),SP,par_io_kind='independent') + call io_variable_bulk(ID_frag,1,C3=MPA_E_par(1)%blc(:,:,:X%mpa_npoles),IPOS=(/1,MPA_E_par(1)%rows(1),MPA_E_par(1)%cols(1),1/)) ! ! residues @iq ! ch="MPA_R_Q_"//trim(intc(iq)) ! - if(write_is_on(ID)) then -#if defined _PAR_IO - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP,par_io_kind='collective') -#else - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP,par_io_kind='independent') -#endif - call io_variable_bulk(ID_frag,1,C3=MPA_R_par(1)%blc(:,:,:Xw%n_freqs),IPOS=(/1,MPA_R_par(1)%rows(1),MPA_R_par(1)%cols(1),1/)) - ! - else if(read_is_on(ID)) then - call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,Xw%n_freqs/),SP) - call io_variable_bulk(ID_frag,1,C3=MPA_R_par(1)%blc(:,:,:Xw%n_freqs),IPOS=(/1,MPA_R_par(1)%rows(1),MPA_R_par(1)%cols(1),1/)) - endif + call def_variable_bulk(ID_frag,trim(ch),1,(/2,X%ng_db,X%ng_db,X%mpa_npoles/),SP,par_io_kind='independent') + call io_variable_bulk(ID_frag,1,C3=MPA_R_par(1)%blc(:,:,:X%mpa_npoles),IPOS=(/1,MPA_R_par(1)%rows(1),MPA_R_par(1)%cols(1),1/)) ! if (read_is_on(ID) .and. different_db_RL_order) then call error('[io_MPA] different_db_RL_order not implemented') diff --git a/src/modules/SET_defaults.F b/src/modules/SET_defaults.F index 97c94566e0..04e2842a1b 100644 --- a/src/modules/SET_defaults.F +++ b/src/modules/SET_defaults.F @@ -19,7 +19,7 @@ subroutine SET_defaults(INSTR,IND,OD,COM_DIR) & l_X_terminator,X_terminator_E,global_gauge,& & Chi_linalg_mode,X_use_lin_sys,X_use_gpu,X_DbGd_percentual use QP_m, ONLY:QP_dSc_steps,QP_n_W_freqs,QP_G_Zoom_treshold,& -& QP_dSc_test,QP_solver,QP_G_damp,QP_dSc_delta,& +& QP_dSc_test,QP_solver,QP_G_damp,QP_G_mprop,QP_dSc_delta,& & QP_cg_percent,QP_n_states,SC_E_threshold, & & QP_Sc_steps,QP_G_er,QP_G_ir,QP_G_dr,SC_band_mixing,QP_G_solver,& & COHSEX_use_empties,On_Mass_Shell_approx,& @@ -285,7 +285,7 @@ subroutine SET_defaults(INSTR,IND,OD,COM_DIR) ! RIM-W ! RIM_W_ng=0 - RIMW_type="default" + rimw_type='semiconductor' ! ! CUTOFF ! @@ -374,6 +374,7 @@ subroutine SET_defaults(INSTR,IND,OD,COM_DIR) QP_dSc_test=.FALSE. QP_solver=' ' QP_G_damp=0.1/HA2EV + QP_G_mprop=1 QP_dSc_delta=0.1/HA2EV QP_G_solver='ra' QP_G_er=(/-10._SP/HA2EV,10._SP/HA2EV/) diff --git a/src/modules/mod_QP.F b/src/modules/mod_QP.F index cb786e2d9f..ee27775811 100644 --- a/src/modules/mod_QP.F +++ b/src/modules/mod_QP.F @@ -42,6 +42,7 @@ module QP_m real(SP) :: QP_G_Zoom_treshold real(SP) :: QP_time_order_sign=-1 ! T-ordered. Set to +1 (causal) in QP_SET_temperature_pre_factor real(SP) :: QP_G_damp + integer :: QP_G_mprop ! GS: order of Lorentzians used in the G[W] propagator real(SP) :: QP_G_er(2) real(SP) :: QP_G_ir(2) real(SP) :: QP_G_dr(2) @@ -65,7 +66,7 @@ module QP_m logical, allocatable:: QP_state(:,:) ! ! Solver... - character(2) :: QP_G_solver + character(2) :: QP_G_solver ! DALV: for G&Sigma_mpa real(SP) :: QP_dSc_delta logical :: QP_dSc_test logical :: On_Mass_Shell_approx diff --git a/src/modules/mod_R_lattice.F b/src/modules/mod_R_lattice.F index 0dafcb238d..64e0d96e13 100644 --- a/src/modules/mod_R_lattice.F +++ b/src/modules/mod_R_lattice.F @@ -153,9 +153,17 @@ module R_lattice logical :: RIM_W_is_diagonal real(SP), allocatable, target :: RIM_W_E(:,:) complex(SP),allocatable, target :: RIM_W (:,:,:,:) - complex(SP),allocatable, target DEV_ATTR :: RIM_W_d(:,:,:,:) - complex(DP),allocatable :: f_coeff(:,:,:,:,:) - character(schlen) :: RIMW_type + complex(SP),allocatable, target DEV_ATTR :: RIM_W_d (:,:,:,:) + complex(SP),allocatable, target :: aux_RIM_W(:,:,:) + ! AF: not needed at present + !complex(SP),allocatable, target DEV_ATTR :: aux_RIM_W_d(:,:,:) + complex(DP),allocatable, target :: f_coeff(:,:,:,:,:) + complex(DP),allocatable, target DEV_ATTR :: f_coeff_d(:,:,:,:,:) + character(schlen) :: rimw_type, rimw_check + ! + ! System-dependent screening terms + logical :: l_Dirac + real(SP) :: Dirac_v_f ! Dirac cone Fermi velocity (x 10^6 m/s) ! ! Coulomb (including Cutoff) ! diff --git a/src/modules/mod_X.F b/src/modules/mod_X.F index 2fd85b9bc8..453516c010 100644 --- a/src/modules/mod_X.F +++ b/src/modules/mod_X.F @@ -124,8 +124,8 @@ module X_m ! type X_t ! DS: is 1 Xo (see X_dielectric_matrix:109) or Xx (see io_X.F) ? - ! why is 5 reported as IP in the develop - integer :: whoami ! 1:Xo/Xx 2:em1s 3:em1d 4:ppa 5:mpa ?:IP + ! + integer :: whoami ! 1:Xo/Xx 2:em1s 3:em1d 4:ppa 5:mpa integer :: ng integer :: ng_db integer :: iq(2) diff --git a/src/modules/mod_memory.F b/src/modules/mod_memory.F index 8852053b70..ba91e40561 100644 --- a/src/modules/mod_memory.F +++ b/src/modules/mod_memory.F @@ -521,8 +521,9 @@ subroutine MEM_rd6_d(name,c) real(DP) DEV_ATTR ::c(:,:,:,:,:,:) call MEM_manager_alloc(name,size(c,KIND=IPL),kind(c(1,1,1,1,1,1)),DEV_) end subroutine - + ! #endif + ! ! end of GPU subroutines ! subroutine MEM_l1(name,l) diff --git a/src/parser/mod_it_tools.F b/src/parser/mod_it_tools.F index 7dd694d633..7a475a4004 100644 --- a/src/parser/mod_it_tools.F +++ b/src/parser/mod_it_tools.F @@ -87,16 +87,17 @@ subroutine it_reset(mode) logical function check_verbosity(what) use drivers, ONLY:infile_editing character(*) :: what - if (what=="kpt") check_verbosity=infile_verbosity==V_kpt.or.infile_verbosity==V_all - if (what=="RL" ) check_verbosity=infile_verbosity==V_RL.or.infile_verbosity==V_all - if (what=="sc" ) check_verbosity=infile_verbosity==V_sc.or.infile_verbosity==V_all - if (what=="qp" ) check_verbosity=infile_verbosity==V_qp.or.infile_verbosity==V_all - if (what=="io" ) check_verbosity=infile_verbosity==V_io.or.infile_verbosity==V_all + if (what=="kpt") check_verbosity=infile_verbosity==V_kpt.or.infile_verbosity==V_all + if (what=="RL" ) check_verbosity=infile_verbosity==V_RL.or.infile_verbosity==V_all + if (what=="sc" ) check_verbosity=infile_verbosity==V_sc.or.infile_verbosity==V_all + if (what=="qp" ) check_verbosity=infile_verbosity==V_qp.or.infile_verbosity==V_all + if (what=="io" ) check_verbosity=infile_verbosity==V_io.or.infile_verbosity==V_all if (what=="general" ) check_verbosity=infile_verbosity==V_general.or.infile_verbosity==V_all if (what=="resp" ) check_verbosity=infile_verbosity==V_resp.or.infile_verbosity==V_all if (what=="real_time") check_verbosity=infile_verbosity==V_real_time.or.infile_verbosity==V_all if (what=="parallel" ) check_verbosity=infile_verbosity==V_parallel.or.infile_verbosity==V_all - if (what=="nl" ) check_verbosity=infile_verbosity==V_nl_optics.or.infile_verbosity==V_all + if (what=="nl" ) check_verbosity=infile_verbosity==V_nl_optics.or.infile_verbosity==V_all + if (what=="coul" ) check_verbosity=infile_verbosity==V_coul.or.infile_verbosity==V_all if (.not.infile_editing) check_verbosity=.TRUE. end function ! diff --git a/src/parser/mod_itm.F b/src/parser/mod_itm.F index a42d011529..d01585fabb 100644 --- a/src/parser/mod_itm.F +++ b/src/parser/mod_itm.F @@ -38,6 +38,7 @@ module it_m integer, parameter :: V_parallel=9 integer, parameter :: V_nl_optics=10 integer, parameter :: V_ph=11 + integer, parameter :: V_coul=12 integer, parameter :: V_all=99 ! integer, parameter :: maxrnlvls=60 diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index eab77a26a5..aa0953268b 100644 --- a/src/pol_function/X_dielectric_matrix.F +++ b/src/pol_function/X_dielectric_matrix.F @@ -20,18 +20,20 @@ integer function X_dielectric_matrix(Xen,Xk,q,X,Xw,Dip,SILENT_MODE,CHILD) use LIVE_t, ONLY:live_timing use stderr, ONLY:intc use frequency, ONLY:w_samp - use R_lattice, ONLY:bz_samp,nqibz + use R_lattice, ONLY:bz_samp,nqibz,l_Dirac,Dirac_v_f use electrons, ONLY:levels use parallel_int, ONLY:PP_wait,PARALLEL_global_indexes,PARALLEL_WF_distribute,PARALLEL_WF_index use parallel_m, ONLY:PAR_IND_Q_ibz,PAR_Q_ibz_index,PAR_nQ_ibz,PAR_COM_X_WORLD,& & PAR_IND_Xk_ibz,PAR_IND_CON_BANDS_X,PAR_IND_VAL_BANDS_X use wave_func, ONLY:WF_buffered_IO,WF,WF_buffer use IO_int, ONLY:io_control,IO_and_Messaging_switch - use IO_m, ONLY:OP_RD_CL,OP_APP_CL,VERIFY,REP,io_RESPONSE + use IO_m, ONLY:OP_RD_CL,OP_APP_CL,VERIFY,REP,io_RESPONSE,OP_WR_CL + use com, ONLY:depth,msg use TDDFT, ONLY:F_xc_gspace use interfaces, ONLY:WF_load,WF_free use QP_m, ONLY:QP_n_W_freqs_redux use matrix, ONLY:MATRIX_reset,MATRIX_copy + use parser_m, ONLY:parser ! #include ! @@ -73,6 +75,15 @@ integer function X_dielectric_matrix(Xen,Xk,q,X,Xw,Dip,SILENT_MODE,CHILD) ! call TDDFT_do_X_W_typs(-1,X,Xw) ! + !Check if to include the Dirac cone Fermi velocity + ! + call parser('graphene',l_Dirac) + if(l_Dirac) then + call msg('r','Graphene q--> 0 term included') + call parser('v_fermi',Dirac_v_f) + end if + ! + ! ! Sectioning ! sec_mode='*' diff --git a/src/pol_function/X_irredux.F b/src/pol_function/X_irredux.F index 60a9173882..4fe9eee647 100644 --- a/src/pol_function/X_irredux.F +++ b/src/pol_function/X_irredux.F @@ -24,14 +24,14 @@ subroutine X_irredux(iq,what,X_par,Xen,Xk,Xw,X,Dip) use drivers, ONLY:l_life use IO_m, ONLY:io_RESPONSE,io_DIP use ALLOC, ONLY:DIPOLE_ALLOC_global - use pars, ONLY:SP,cZERO,schlen + use pars, ONLY:SP,cZERO,schlen,pi use wrapper, ONLY:V_plus_alpha_V,vv_caxpy use LIVE_t, ONLY:live_timing use com, ONLY:msg use matrix, ONLY:PAR_matrix use stderr, ONLY:intc use wave_func, ONLY:WF - use parallel_m, ONLY:PAR_COM_X_WORLD,PAR_COM_X_WORLD_RL_resolved,myid + use parallel_m, ONLY:PAR_COM_X_WORLD,PAR_COM_RL_INDEX,myid,PAR_COM_X_WORLD_RL_resolved,master_cpu use parallel_int, ONLY:PP_redux_wait use openmp, ONLY:OPENMP_update,n_threads_X,master_thread,OPENMP_set_threads,n_threads_X,& & n_out_threads,n_inn_threads,OPENMP_locks_reset,n_threads_FFT @@ -41,8 +41,8 @@ subroutine X_irredux(iq,what,X_par,Xen,Xk,Xw,X,Dip) use frequency, ONLY:w_samp,bare_grid_N,coarse_grid_N,coarse_grid_Pt use interfaces, ONLY:WF_load use electrons, ONLY:levels - use R_lattice, ONLY:qindx_X,bz_samp,G_m_G - use D_lattice, ONLY:i_space_inv + use R_lattice, ONLY:qindx_X,bz_samp,G_m_G,q_norm,idir,l_Dirac,Dirac_v_f + use D_lattice, ONLY:i_space_inv,alat use collision_el, ONLY:elemental_collision,elemental_collision_free,elemental_collision_alloc use DIPOLES, ONLY:DIPOLE_t use X_m, ONLY:X_t,X_poles,X_Ein_poles,current_iq,X_poles_tab,& @@ -55,6 +55,7 @@ subroutine X_irredux(iq,what,X_par,Xen,Xk,Xw,X,Dip) use gpu_m, ONLY:have_gpu use devxlib, ONLY:devxlib_memcpy_d2d,devxlib_memcpy_d2h,devxlib_memcpy_h2d use timing_m, ONLY:timing + use units, ONLY:HA2EV ! #include #include @@ -74,8 +75,8 @@ subroutine X_irredux(iq,what,X_par,Xen,Xk,Xw,X,Dip) integer :: ig1,ig2,iw,n_poles,i_cg,i_bg,mutexid,ngrho,& & X_rows1,X_rows2,X_cols1,X_cols2,X_nrows,io_err logical :: force_bare_X_G,Drude_pole,skip_WF_load - real(SP) :: minmax_ehe(2,PAR_COM_X_WORLD%n_CPU),cutoff - complex(SP) :: GreenF(Xw%n_freqs),drude_GreenF(Xw%n_freqs),ctmp(1,1) + real(SP) :: minmax_ehe(2,PAR_COM_X_WORLD%n_CPU),cutoff,lcut,Dirac_v_f_au + complex(SP) :: GreenF(Xw%n_freqs),drude_GreenF(Xw%n_freqs),ctmp(1,1),X_Dirac complex(SP),allocatable, target DEV_ATTR :: Xo_res(:,:) complex(SP),allocatable DEV_ATTR :: work(:) @@ -283,6 +284,40 @@ subroutine X_irredux(iq,what,X_par,Xen,Xk,Xw,X,Dip) enddo endif ! + !Graphene intraband term + ! + if (iq==1.and.master_thread.and.(l_Dirac)) then + ! + if(X_par%rows(1)==1.and.X_par%cols(1)==1) then + !call msg('r', 'Graphene analytical model of intraband transitions') + call warning('Dirac cone contribution is applied only to head static irr. polarizability') + ! + ! Head correction + !================= + ! + lcut = alat(idir(1)) + ! + call msg('r', 'Fermi velocity (10^6 m/s)', Dirac_v_f) + Dirac_v_f_au = Dirac_v_f/2.187691_SP + call msg('r', 'Fermi velocity (a.u.)', Dirac_v_f_au) + ! + !Analytical model graphene: + !X^0_00_Dirac(w=0) = -1/(4*v_fermi)*q/L + ! + X_Dirac = -q_norm(1)/(4._SP*Dirac_v_f_au*lcut) + ! + if(master_cpu) then + ! + X_par%blc(1,1,1)=X_par%blc(1,1,1)+X_Dirac + ! + if (have_gpu) X_par%blc_d(1,1,1)=X_par%blc(1,1,1) + ! + endif + ! + endif + endif + ! + ! ! MAIN LOOP !=========== ! diff --git a/src/qp/QP_descriptions.F b/src/qp/QP_descriptions.F index fd337695e4..51c6e3c81d 100644 --- a/src/qp/QP_descriptions.F +++ b/src/qp/QP_descriptions.F @@ -12,7 +12,7 @@ subroutine QP_descriptions(qp,X,Xw,Update) use units, ONLY:HA2EV use X_m, ONLY:X_t,use_X_DbGd,OPTICAL_averaged_dirs_string use QP_m, ONLY:QP_t,QP_ng_Sx,QP_n_G_bands,QP_cg_percent,& -& QP_dSc_delta,QP_G_damp,QP_dSc_steps,& +& QP_dSc_delta,QP_G_damp,QP_G_mprop,QP_dSc_steps,& & COHSEX_use_empties,QP_Sc_steps,QP_G_er,QP_G_dr,QP_solver,& & use_GreenF_to_eval_QP,GF_is_causal,QP_G_Zoom_treshold,& & GW_terminator_E,l_GW_terminator,GW_terminator_Kind @@ -121,6 +121,7 @@ subroutine QP_descriptions(qp,X,Xw,Update) endif else call IO_desc_add(qp%desc,str='Sc/G damping',term='ev',kind='r',R=(/QP_G_damp*HA2EV/)) + call IO_desc_add(qp%desc,str='Sc/G propagator order',kind='i',I=(/QP_G_mprop/)) endif ! call IO_desc_add(qp%desc,str='Sc bands terminator',kind='l',L=l_GW_terminator) diff --git a/src/qp/QP_driver.F b/src/qp/QP_driver.F index 9cd663fd19..588f7a2e33 100644 --- a/src/qp/QP_driver.F +++ b/src/qp/QP_driver.F @@ -305,8 +305,15 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw,Dip) ! ! GF for ppa and mpa ! - if ( l_ppa ) call QP_ppa_cohsex(X,Xk,en,k,q,qp,Xw,0) - if ( l_mpa ) call QP_mpa(X,Xk,en,k,q,qp,Xw,0) + if ( l_ppa ) then + if (l_rim_w) call QP_interpolate_W(X,Xw,q,'PPA') + call QP_ppa_cohsex(X,Xk,en,k,q,qp,Xw,0) + end if + ! + if ( l_mpa ) then + if (l_rim_w) call QP_interpolate_W(X,Xw,q,'MPA') + call QP_mpa(X,Xk,en,k,q,qp,Xw,0) + end if ! endif ! diff --git a/src/qp/QP_interpolate_W.F b/src/qp/QP_interpolate_W.F index 78a79b55c1..2792b0e2bf 100644 --- a/src/qp/QP_interpolate_W.F +++ b/src/qp/QP_interpolate_W.F @@ -3,15 +3,20 @@ ! ! Copyright (C) 2006 The Yambo Team ! -! Authors (see AUTHORS file for details): AG +! Authors (see AUTHORS file for details): AG GS +! +! AF: This source needs cleanup +! - proper indentation +! - del of old/unused source code +! - ... ! subroutine QP_interpolate_W(X,Xw,q,mode) ! use pars, ONLY:SP,pi,zero_dfl,schlen,DP,rZERO use com, ONLY:msg - use R_lattice, ONLY:bz_samp,RIM_W_ng,b,k_grid_b,bare_qpg,& -& RIM_W_is_diagonal,RIM_W,RIM_W_d,f_coeff,idir,RIM_W_E,& -& RIM_id_epsm1_reference,RIM_epsm1,RIM_qpg + use R_lattice, ONLY:bz_samp,RIM_W_ng,g_vec,b,k_grid_b,bare_qpg,& +& RIM_W_is_diagonal,RIM_W,f_coeff,idir,RIM_W_E,cut_is_slab,& +& RIM_id_epsm1_reference,RIM_epsm1,Dirac_v_f,RIM_qpg,rimW_type use vec_operate, ONLY:c2a,v_norm use ALLOC, ONLY:X_ALLOC_elemental use X_m, ONLY:X_mat,X_t @@ -23,6 +28,7 @@ subroutine QP_interpolate_W(X,Xw,q,mode) use D_lattice, ONLY:alat,a use timing_m, ONLY:timing use parallel_m, ONLY:master_cpu + use parser_m, ONLY:parser ! #include #include @@ -34,54 +40,74 @@ subroutine QP_interpolate_W(X,Xw,q,mode) ! ! Work Space ! - integer :: iq1,ig1,iq2,ig2,iq,ig + integer :: iq1,ig1,iq2,ig2,iq,ig,i1,i2,iq0,iq3,iqn1,iqn2 + integer :: ig1_tmp,iq2_tmp,ig2_tmp integer :: nn,idm,igr,igc,ig2max,iomega - integer :: idx_q(RIM_W_ng,q%nibz,5),idx_G(RIM_W_ng,q%nibz,5),idx_is(RIM_W_ng,q%nibz,5) + integer :: idx_q(RIM_W_ng,q%nibz,7),idx_G(RIM_W_ng,q%nibz,7),idx_is(RIM_W_ng,q%nibz,7) + integer :: idx_q_2(RIM_W_ng,q%nibz,13),idx_G_2(RIM_W_ng,q%nibz,13),idx_is_2(RIM_W_ng,q%nibz,13) real(SP) :: dummy(2),anis_factor,em1_anis(3) - complex(DP) :: f_func(5),dp_dummy(5) - complex(DP) :: vX_nn(5,q%nibz,RIM_W_ng,RIM_W_ng,Xw%n_freqs) + complex(SP) :: f_func(13),dp_dummy(13),dp_temp(3) + complex(DP) :: vX_nn(7,q%nibz,RIM_W_ng,RIM_W_ng,Xw%n_freqs),vX_nnn(13,q%nibz,RIM_W_ng,RIM_W_ng,Xw%n_freqs) real(SP) :: v1(3),r_dum(2) - real(SP) :: q_grid_b_rlu(3,2),q_grid_b_iku(3,2),q_grid_b_cc(3,2) - integer :: ID, IO_ACT, io_err, ng_save,G_max_ibz(q%nibz), G_max + real(SP) :: Dirac_v_f_au,X0_Dirac + real(SP) :: q_grid_b_rlu(3,3),q_grid_b_iku(3,3),q_grid_b_cc(3,3) + integer :: ID, IO_ACT,io_err, ng_save,G_max_ibz(q%nibz), G_max integer, external :: io_X, io_RIM_W character(schlen) :: ch_dum,fmt_dum,msg_dum(2) + logical :: l_check,q0_neighbour,q0_snd_neighbour ! call timing('RIM-W-coeff',OPR='start') ! - YAMBO_ALLOC(f_coeff,(6,RIM_W_ng,RIM_W_ng,q%nibz,Xw%n_freqs)) + YAMBO_ALLOC(f_coeff,(20,RIM_W_ng,RIM_W_ng,q%nibz,Xw%n_freqs)) ! call section("+", "RIM-W interpolation") ! - !Load q_grid_b in iku - do idm=1,2 + call msg("r", trim(rimW_type)) + ! + !Load q_grid_b in rlu + do idm=1,3 v1 = k_grid_b(idm,:) call c2a(v_in=v1,mode='kc2a') ! - if (abs(v1(1)) < zero_dfl) then + if((abs(v1(2)) X%ng) call error(' G_max '//trim(intc(G_max))//& -& ' for RIM-W is higher than G_max '//trim(intc(X%ng))//' of vX') + & ' for RIM-W is higher than G_max '//trim(intc(X%ng))//' of vX') ! call msg('r', 'G-vectors loaded', maxval(G_max_ibz)) call msg('r', 'Number of interpolated frequencies', Xw%n_freqs) - call msg('r','') ! - X%ng = G_max_ibz(1) + ! + !Load only the G needed + X%ng = G_max call io_control(ACTION=OP_RD,COM=REP,SEC=(/1/),ID=ID) io_err=io_X(X,Xw,ID) ! @@ -138,16 +175,30 @@ subroutine QP_interpolate_W(X,Xw,q,mode) ! ! I/O of X ! - !Load only the G needed - X%ng = G_max_ibz(iq) + !call msg("rs","iq =", iq) + !call msg("rs","Gmax", G_max_ibz(iq)) ! IO_ACT=manage_action(RD_CL_IF_END,iq,1,q%nibz) call io_control(ACTION=IO_ACT,COM=NONE,SEC=(/2*iq,2*iq+1/),ID=ID) + !if (iq==1) then + ! call io_control(ACTION=OP_RD,COM=REP, SEC=(/1,2,2*iq+1/),ID=ID) + ! else if (q%nibz==1) then + ! call io_control(ACTION=OP_RD_CL,COM=NONE, SEC=(/1,2,3/),ID=ID) + ! else if (iq > 1) then + ! call io_control(ACTION=RD_CL_IF_END,COM=NONE,SEC=(/2*iq,2*iq+1/),ID=ID) + !endif + ! ! io_err=io_X(X,Xw,ID) ! + ! + !call msg('r', 'Re frequency 1', real(Xw%p(1))) + !call msg('r', 'Im frequency 1', aimag(Xw%p(1))) + ! call deliver_IO_error_message(io_err,'PP/Em1s',STOP_it=.TRUE.) ! + !call msg("rs","shape X_mat",shape(X_mat)) + ! !Assign X do iomega=1,Xw%n_freqs do ig1=1,RIM_W_ng @@ -155,26 +206,80 @@ subroutine QP_interpolate_W(X,Xw,q,mode) ig2max = RIM_W_ng if (RIM_W_is_diagonal) ig2max=ig1 do iq1=1,q%nibz - do ig2=ig1,ig2max +! if (v_norm(q%pt(iq1,:))<0.2) then + do ig2=1,ig2max +! old version do ig2=ig1,ig2max ! - do nn=1,5 + do nn=1,7 ! if (idx_q(ig1,iq1,nn) == iq) then vX_nn(nn,iq1,ig1,ig2,iomega) = X_mat(idx_G(ig1,iq1,nn),idx_G(ig2,iq1,nn),iomega) endif + if (idx_q_2(ig1,iq1,nn) == iq) then + vX_nnn(nn,iq1,ig1,ig2,iomega) = X_mat(idx_G_2(ig1,iq1,nn),idx_G_2(ig2,iq1,nn),iomega) + end if ! enddo - ! - enddo - enddo + ! + ! Calculation of next nearest neighbours vX needed only for the q=0 term and its neighbours + ! in the head term + if (ig1==1 .and. ig2==1) then + do iqn2=1,7 + if (idx_q(1,iq,iqn2)==1) then + do nn=8,13 + if (idx_q_2(ig1,iq1,nn) == iq) then + vX_nnn(nn,iq1,ig1,ig2,iomega) = X_mat(idx_G_2(ig1,iq1,nn),idx_G_2(ig2,iq1,nn),iomega) + end if + end do + cycle + end if + end do + end if + end do + end do + ! + ! Calculation of next nearest neighbours vX needed only for the q=0 term and the neighboursing q-points + ! for the wings and the head term + ! + ! do iqn2=1,7 + ! ! +! if (idx_q(1,iq1,iqn2)==1) then +! do nn=1,13 + ! if (idx_q_2(1,iq1,nn) == iq) then + ! vX_nnn(nn,iq1,1,ig1,iomega) = X_mat(idx_G_2(1,iq1,nn),idx_G_2(ig1,iq1,nn),iomega) + ! end if +! end do + ! end if +! end do + ! +! end if +! end do + ! enddo enddo ! - call X_ALLOC_elemental('X') + call X_ALLOC_global('X') + !if(iq==15) then + ! write(*,*) "Re[X_mat]" + ! write(*,*) real(X_mat(:,:,1)) + ! write(*,*) "Im[X_mat]" + ! write(*,*) aimag(X_mat(:,:,1)) + !end if ! enddo + !call msg("rs","Re[vX] incrimined", real(vX_nn(1,16,1,1,1))) + !call msg("rs","Im[vX] incrimined", aimag(vX_nn(1,16,1,1,1))) + !call msg("rs","Re[vX] incrimined", real(vX_nn(2,16,1,1,1))) + !call msg("rs","Im[vX] incrimined", aimag(vX_nn(2,16,1,1,1))) + !call msg("rs","Re[vX] incrimined", real(vX_nn(3,16,1,1,1))) + !call msg("rs","Im[vX] incrimined", aimag(vX_nn(3,16,1,1,1))) + !call msg("rs","Re[vX] incrimined", real(vX_nn(4,16,1,1,1))) + !call msg("rs","Im[vX] incrimined", aimag(vX_nn(4,16,1,1,1))) + !call msg("rs","Re[vX] incrimined", real(vX_nn(5,16,1,1,1))) + !call msg("rs","Im[vX] incrimined", aimag(vX_nn(5,16,1,1,1))) ! X%ng = ng_save + YAMBO_FREE(X_mat) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !Calculation of the interpolation coefficients @@ -187,87 +292,791 @@ subroutine QP_interpolate_W(X,Xw,q,mode) if (RIM_id_epsm1_reference==0) em1_anis=0. ! do iomega=1,Xw%n_freqs - do igr=1,RIM_W_ng - ! - ig2max = RIM_W_ng - if (RIM_W_is_diagonal) ig2max=igr - do igc=igr,ig2max + !call msg('r','iomega',iomega) + do igr=1,RIM_W_ng ! - do iq1=1,q%nibz + ig2max = RIM_W_ng + if (RIM_W_is_diagonal) ig2max=igr + do igc=1,ig2max + ! old version do igc=igr,ig2max ! - ! Evaluate f_func nearest neighbour - ! - do nn=1,5 - ! - !Select the index of the reference/n.n. - iq2 = idx_q(igr,iq1,nn) - ig1 = idx_G(igr,iq1,nn) - ig2 = idx_G(igc,iq1,nn) - ! - f_func(nn) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& -& vX_nn(nn,iq1,igr,igc,iomega)/(vX_nn(nn,iq1,igr,igc,iomega)+1._DP) - ! - enddo - ! - ! q == 0 terms must be tratened separately - if (iq1 == 1 .and. igr == 1 .and. igc == 1) then - ! - call msg("r","Evaluating coefficients q = 0") + do iq1=1,q%nibz +! if (v_norm(q%pt(iq1,:))<0.2) then + q0_neighbour=.false. + q0_snd_neighbour=.false. + !cut_is_slab + ! Evaluate f_func nearest neighbour ! - if (RIM_id_epsm1_reference/=0 .and. iomega==1) then + do nn=1,7 ! - !vX is included in the anis_factor - call msg('r','This is the anysotropy case') - f_coeff(1,igr,igc,iq1,iomega)=(4._DP*pi)/(bare_qpg(iq1,igr)*bare_qpg(iq1,igc)*(2._DP*pi*alat(idir(1)))**2) + !Select the index of the reference/n.n. + iq2 = idx_q(igr,iq1,nn) + ig1 = idx_G(igr,iq1,nn) + ig2 = idx_G(igc,iq1,nn) ! - else - f_coeff(1,igr,igc,iq1,iomega) = vX_nn(1,iq1,igr,igc,iomega)*(4._DP*pi) / & -& (bare_qpg(iq1,igr)*bare_qpg(iq1,igc)*(2._DP*pi*alat(idir(1)))**2) - end if - ! - do nn=2,3 - idm = MOD(nn+1,2)+1 + f_func(nn) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nn(nn,iq1,igr,igc,iomega)/(vX_nn(nn,iq1,igr,igc,iomega)+1._DP) + ! + ! For metals, special treatment if a neighbour is the q-->0 term at zero frequency, + ! as this term is ill defined in metals both in the case of the head and wings + ! + if ((iq1==1) .and. (iq2==1) .and. (ig1==1) .and. (ig2 .ne. 1) .and. (iomega .ne. 1)) then + f_func(nn)=0 + end if ! - if (RIM_id_epsm1_reference/=0 .and. iomega==1) then + ! The q-->0 term is a neighbour + if ( (iq1 .ne. 1) .and. (iq2==1)) then + select case (trim(rimw_type)) + case ("semiconductor") + ! + if ( (ig1==1) .and. (ig2==1)) then + if (.not. cut_is_slab) then + q0_neighbour=.true. + iq0=nn + f_func(nn)=vX_nn(1,1,1,1,iomega)/(4._DP*pi*(1._DP+vX_nn(1,1,1,1,iomega))) + end if + end if ! - anis_factor=dot_product(em1_anis,q_grid_b_cc(:,idm)**2) - f_coeff(nn,igr,igc,iq1,iomega)=-1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & -& (f_coeff(1,iq1,igc,iq1,iomega)*anis_factor)) - else - f_coeff(nn,igr,igc,iq1,iomega)=-1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & -& (v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,iq1,igc,iq1,iomega))) + case ("metal") + ! Head + if ( (ig1==1) .and. (ig2==1)) then + if (iomega==1) then + f_func(nn)=f_coeff(1,1,1,1,1) + else + f_func(nn)=0 + q0_neighbour=.true. + iq0=nn + end if + ! Wing + else if ((ig1==1) .and. (ig2 .ne. 1) ) then + f_func(nn)=f_coeff(1,ig1,ig2,1,iomega) + if (iomega==1) then + f_func(nn)=0 + if (.not. cut_is_slab) then + q0_neighbour=.true. + iq0=nn + end if + end if + end if + end select end if ! + ! + if ((ig1==1) .and. (ig2==1) .and. (iomega .ne. 1) .and. (iq2 .ne. 1) .and. (iq1 .ne. 1)) then + iq3 = idx_q(ig1,iq2,nn) + if (iq3==1) then + q0_snd_neighbour=.true. + iq0=nn + end if + end if enddo ! - else + !test = (iq1==1.and.igr == 1.and.igc == 1) + !call msg('r','iq',iq1) + !call msg('r','igr',igr) + !call msg('r','igc',igc) + !call msg('r','test',test) ! - !call msg("r","Evaluating coefficients q/= 0") + ! Coefficients calculated from Taylor expansion + ! Coefficients in rlu ! - !Coefficients in rlu + dp_dummy(:) = 0 dp_dummy(1) = f_func(1) !f0 - dp_dummy(2) = (f_func(3)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)) !f1 - dp_dummy(3) = (f_func(2)-f_func(4))/(2._DP*q_grid_b_rlu(2,2)) !f2 - dp_dummy(4) = (f_func(3)-2._DP*f_func(1)+f_func(5))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 - dp_dummy(5) = (f_func(2)-2._DP*f_func(1)+f_func(4))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + dp_dummy(2) = (f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)) !f1 + dp_dummy(3) = (f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)) !f2 + dp_dummy(4) = (f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)) !f3 + dp_dummy(5) = (f_func(2)-2._DP*f_func(1)+f_func(5))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 + dp_dummy(6) = (f_func(3)-2._DP*f_func(1)+f_func(6))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + dp_dummy(7) = (f_func(4)-2._DP*f_func(1)+f_func(7))/(2._DP*q_grid_b_rlu(3,3)**2) !f33 + ! + ! for 3D + ! if nearest neighbour is q=0 along a certain direction, use cubic expansion of the f functions + ! when the lowest order is quadratic + ! + if (q0_neighbour) then + dp_dummy(2:4)=0 + if (iq0>=5) then + i1=iq0-3 + dp_temp(1)= 0 + dp_temp(2)= (8*f_func(1)-f_func(i1))/(4*q_grid_b_rlu(i1-1,i1-1)**2) + dp_temp(3)= (f_func(i1)-4*f_func(1))/(4*q_grid_b_rlu(i1-1,i1-1)**3) + ! + dp_dummy(1) = dp_temp(1)*q_grid_b_rlu(i1-1,i1-1)+dp_temp(2)*q_grid_b_rlu(i1-1,i1-1)**2 + & +& dp_temp(3)*q_grid_b_rlu(i1-1,i1-1)**3 + dp_dummy(i1) = dp_temp(1)+2*dp_temp(2)*q_grid_b_rlu(i1-1,i1-1) + & +& 3*dp_temp(3)*q_grid_b_rlu(i1-1,i1-1)**2 + dp_dummy(i1+3) = dp_temp(2)+3*dp_temp(3)*q_grid_b_rlu(i1-1,i1-1) + dp_dummy(i1+6) = dp_temp(3) + ! + ! a mixing quadratic term qi1*qi2 is needed in case the i1 direction is not orthogonal to the other directions + do idm=i1-1,i1 + dp_dummy(8+i1+mod(idm,3)) = (f_func(2+mod(idm,3))-f_func(5+mod(idm,3)))/ & +& (4._DP*(q_grid_b_rlu(i1-1,i1-1)*q_grid_b_rlu(1+mod(idm,3),1+mod(idm,3)))) + dp_dummy(2+mod(idm,3)) = 2.0*dp_dummy(8+i1+mod(idm,3))*q_grid_b_rlu(i1-1,i1-1) + end do + ! + ! in case the i1 direction orthogonal to the other two directions, but the direction i2 and i3 are not between + ! each other a mixing quadratic term of the kind qi2*qi3 is needed + ! + if ( (q_grid_b_rlu(i1-1,1+mod(i1-1,3))==0) .and. (q_grid_b_rlu(1+mod(i1-1,3),i1-1)==0) .and. & +& (q_grid_b_rlu(i1-1,1+mod(i1,3))==0) .and. (q_grid_b_rlu(1+mod(i1,3),i1-1)==0) .and. & +& (q_grid_b_rlu(1+mod(i1-1,3),1+mod(i1,3)) .ne. 0) .or. (q_grid_b_rlu(1+mod(i1,3),1+mod(i1-1,3)) .ne. 0)) then + ! + iqn1 = idx_q_2(1,iq1,12-i1) + ! + f_func(5+mod(i1-1,3))= real(bare_qpg(iqn1,1)*bare_qpg(iqn1,1),kind=DP)/(4._DP*pi)* & +& vX_nnn(12-i1,iq1,1,1,iomega)/(vX_nnn(12-i1,iq1,1,1,iomega)+1._DP) + ! + iqn2 = idx_q_2(1,iq1,15-i1) + ! + f_func(5+mod(i1,3))= real(bare_qpg(iqn2,1)*bare_qpg(iqn2,1),kind=DP)/(4._DP*pi)* & +& vX_nnn(15-i1,iq1,1,1,iomega)/(vX_nnn(15-i1,iq1,1,1,iomega)+1._DP) + ! + dp_dummy(15-i1) = (f_func(5+mod(i1-1,3))-f_func(5+mod(i1,3)))/ & +& (4._DP*(q_grid_b_rlu(1+mod(i1,3),1+mod(i1,3))*q_grid_b_rlu(1+mod(i1-1,3),1+mod(i1-1,3)))) + end if + else + i1=iq0 + dp_temp(1)=0 + dp_temp(2)= (8*f_func(1)-f_func(iq0+3))/(4*q_grid_b_rlu(i1-1,i1-1)**2) + dp_temp(3)= -(f_func(iq0+3)-4*f_func(1))/(4*q_grid_b_rlu(i1-1,i1-1)**3) + dp_dummy(1) = -dp_temp(1)*q_grid_b_rlu(i1-1,i1-1)+dp_temp(2)*q_grid_b_rlu(i1-1,i1-1)**2 - & +& dp_temp(3)*q_grid_b_rlu(i1-1,i1-1)**3 + dp_dummy(i1) = dp_temp(1)-2*dp_temp(2)*q_grid_b_rlu(i1-1,i1-1) + & +& 3*dp_temp(3)*q_grid_b_rlu(i1-1,i1-1)**2 + dp_dummy(i1+3) = dp_temp(2)-3*dp_temp(3)*q_grid_b_rlu(i1-1,i1-1) + dp_dummy(i1+6) = dp_temp(3) + ! + ! in case the i1 direction is not orthogonal to the other two directions, a mixing term is needed + do idm=i1-1,i1 + dp_dummy(8+i1+mod(idm,3)) = (f_func(2+mod(idm,3))-f_func(5+mod(idm,3)))/ & +& (4._DP*(q_grid_b_rlu(i1-1,i1-1)*q_grid_b_rlu(1+mod(idm,3),1+mod(idm,3)))) + dp_dummy(2+mod(idm,3)) = -2.0*dp_dummy(8+i1+mod(idm,3))*q_grid_b_rlu(i1-1,i1-1) + end do + ! + ! in case the i1 direction orthogonal to the other two directions, but the direction i2 and i3 are not between + ! each other a mixing quadratic term of the kind qi2*qi3 is needed + ! + if ( (q_grid_b_rlu(i1-1,1+mod(i1-1,3))==0) .and. (q_grid_b_rlu(1+mod(i1-1,3),i1-1)==0) .and. & +& (q_grid_b_rlu(i1-1,1+mod(i1,3))==0) .and. (q_grid_b_rlu(1+mod(i1,3),i1-1)==0) .and. & +& (q_grid_b_rlu(1+mod(i1-1,3),1+mod(i1,3)) .ne. 0) .or. (q_grid_b_rlu(1+mod(i1,3),1+mod(i1-1,3)) .ne. 0)) then + ! + iqn1 = idx_q_2(1,iq1,12-i1) + ! + f_func(5+mod(i1-1,3))= real(bare_qpg(iqn1,1)*bare_qpg(iqn1,1),kind=DP)/(4._DP*pi)* & +& vX_nnn(12-i1,iq1,1,1,iomega)/(vX_nnn(12-i1,iq1,1,1,iomega)+1._DP) + ! + iqn2 = idx_q_2(1,iq1,15-i1) + ! + f_func(5+mod(i1,3))= real(bare_qpg(iqn2,1)*bare_qpg(iqn2,1),kind=DP)/(4._DP*pi)* & +& vX_nnn(15-i1,iq1,1,1,iomega)/(vX_nnn(15-i1,iq1,1,1,iomega)+1._DP) + ! + dp_dummy(15-i1) = (f_func(5+mod(i1-1,3))-f_func(5+mod(i1,3)))/ & +& (4._DP*(q_grid_b_rlu(1+mod(i1,3),1+mod(i1,3))*q_grid_b_rlu(1+mod(i1-1,3),1+mod(i1-1,3)))) + end if + end if + end if + ! + if (q0_snd_neighbour) then + dp_dummy(2:4)=0 + dp_dummy(5) = (f_func(2)-2._DP*f_func(1)+f_func(5))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 + dp_dummy(6) = (f_func(3)-2._DP*f_func(1)+f_func(6))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + dp_dummy(7) = (f_func(4)-2._DP*f_func(1)+f_func(7))/(2._DP*q_grid_b_rlu(3,3)**2) !f33 + dp_dummy(8) = (f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)**3) !f111 + dp_dummy(9) = (f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)**3) !f222 + dp_dummy(10) = (f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)**3) !f333 + if (iq0 >=5) then + i1=iq0-3 + dp_dummy(i1) = (-6*f_func(iq0)+3._DP*f_func(1)+2*f_func(i1))/(6._DP*q_grid_b_rlu(i1-1,i1-1)) + dp_dummy(i1+3)= (f_func(iq0)-2._DP*f_func(1)+f_func(i1))/(2._DP*q_grid_b_rlu(i1-1,i1-1)**2) + dp_dummy(i1+6)= (3*f_func(iq0)-3._DP*f_func(1)+f_func(i1))/(6._DP*q_grid_b_rlu(i1-1,i1-1)**3) + else + i1=iq0 + dp_dummy(i1) = (-6*f_func(iq0)+3._DP*f_func(1)+2*f_func(i1+3))/(6._DP*q_grid_b_rlu(i1-1,i1-1)) + dp_dummy(i1+3)= (f_func(iq0)-2._DP*f_func(1)+f_func(i1+3))/(2._DP*q_grid_b_rlu(i1-1,i1-1)**2) + dp_dummy(i1+6)= (3*f_func(iq0)-3._DP*f_func(1)+f_func(i1+3))/(6._DP*q_grid_b_rlu(i1-1,i1-1)**3) + end if + end if + ! + if (iq1==1) then + ! + ! Head + !====== + ! + + if ( (igr==1.and.igc==1) & +& .or. ( (igr==1) .and. (igc/= 1) .and. (cut_is_slab) & +& .and. (g_vec(igc,idir(2))==g_vec(igc,idir(3))) .and. (g_vec(igc,idir(2))==0) .and. (iomega==1) ) & +& .or. ( (igr==1) .and. (igc .ne. 1) .and. (cut_is_slab) .and. (g_vec(igc,idir(2))==g_vec(igc,idir(3))) & +& .and. (g_vec(igc,idir(2))==0) .and. (iomega==1) ) ) then + ! + ! Special treatment of the q-->0 term + !===================================== + ! Evaluation of f_func of q-->0 term uses next nearest neighbour + ! Index ordering of NEXT nearest-neighbour maps + ! + ! 6 7 + ! ^ / + ! | / + ! 3 4 + ! ^ / + ! |/ + ! 1 --> 2 --> 5 + ! + ! + idx_q_2 (igc,iq1,1) = iq1 + idx_is_2(igc,iq1,1) = 1 + idx_G_2 (igc,iq1,1) = igc + ! + !Calculate qpG and qpG nearest neighbour + !qpG(:,1) = q%pt(iq1,:)+g_vec(ig1,:) + ! + ! Find nearest 1 --> 2 + ! + iq2_tmp = idx_q(igr,iq1,2) + ig1_tmp = idx_G(igr,iq1,2) + ig2_tmp = idx_G(igc,iq1,2) + ! + ! Find next nearest 2 --> 5 + ! + iq2 = idx_q_2(igr,iq1,5) + ig1 = idx_G_2(igr,iq1,5) + ig2 = idx_G_2(igc,iq1,5) + ! + !call msg('r','iq =',iq1) + !call msg('r','igr =',igr) + !call msg('r','igc =',igc) + ! + !call msg('r','iq_nn =',iq2_tmp) + !call msg('r','igr_nn =',ig1_tmp) + !call msg('r','igc_nn =',ig2_tmp) + !call msg('r','q_nn =',q%pt(iq2_tmp,:)) + ! + !call msg('r','iq_nnn =',iq2) + !call msg('r','igr_nnn =',ig1) + !call msg('r','igc_nnn =',ig2) + !call msg('r','q_nnn =',q%pt(iq2,:)) + ! + f_func(5) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nnn(5,iq1,igr,igc,iomega)/(vX_nnn(5,iq1,igr,igc,iomega)+1._DP) + ! + !Find nearest 1 --> 3 + ! + iq2_tmp = idx_q(igr,iq1,3) + ig1_tmp = idx_G(igr,iq1,3) + ig2_tmp = idx_G(igc,iq1,3) + ! + !Find next nearest 3 --> 6 + ! + iq2 = idx_q_2(igr,iq1,6) + ig1 = idx_G_2(igr,iq1,6) + ig2 = idx_G_2(igc,iq1,6) + ! + !call msg('r','iq_nn =',iq2_tmp) + !call msg('r','igr_nn =',ig1_tmp) + !call msg('r','igc_nn =',ig2_tmp) + !call msg('r','q_nn =',q%pt(iq2_tmp,:)) + ! + !call msg('r','iq_nnn =',iq2) + !call msg('r','igr_nnn =',ig1) + !call msg('r','igc_nnn =',ig2) + !call msg('r','q_nnn =',q%pt(iq2,:)) + ! + f_func(6) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nnn(6,iq1,igr,igc,iomega)/(vX_nnn(6,iq1,igr,igc,iomega)+1._DP) + ! + !Find nearest 1 --> 4 + ! + iq2_tmp = idx_q(igr,iq1,4) + ig1_tmp = idx_G(igr,iq1,4) + ig2_tmp = idx_G(igc,iq1,4) + ! + !Find next nearest 4 --> 7 + ! + iq2 = idx_q_2(igr,iq1,7) + ig1 = idx_G_2(igr,iq1,7) + ig2 = idx_G_2(igc,iq1,7) + ! + !call msg('r','iq_nn =',iq2_tmp) + !call msg('r','igr_nn =',ig1_tmp) + !call msg('r','igc_nn =',ig2_tmp) + !call msg('r','q_nn =',q%pt(iq2_tmp,:)) + ! + !call msg('r','iq_nnn =',iq2) + !call msg('r','igr_nnn =',ig1) + !call msg('r','igc_nnn =',ig2) + !call msg('r','q_nnn =',q%pt(iq2,:)) + ! + f_func(7) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nnn(7,iq1,igr,igc,iomega)/(vX_nnn(7,iq1,igr,igc,iomega)+1._DP) + end if + if (igr==1.and.igc==1) then + if(iomega==1) then + select case (trim(rimw_type)) + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("semiconductor") + ! + call msg("r","Semiconducting q --> 0 shape of f") + ! + if (cut_is_slab) then + !Save f as VXV + if(RIM_id_epsm1_reference/=0) then + ! + !vX is included in the anis_factor + ! + f_coeff(1,1,1,1,1)=(4._DP*pi)/(bare_qpg(1,1)*bare_qpg(1,1)*(2._DP*pi*alat(idir(1)))**2) + else + f_coeff(1,1,1,1,1)=vX_nn(1,1,1,1,1)*(4._DP*pi) / & +& (bare_qpg(1,1)*bare_qpg(1,1)*(2._DP*pi*alat(idir(1)))**2) + end if + ! + do nn=2,4 + idm = nn-1 + ! + if(RIM_id_epsm1_reference/=0 .and. iomega==1) then + anis_factor=dot_product(em1_anis,q_grid_b_cc(:,idm)**2) + f_coeff(nn,1,1,1,1)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (f_coeff(1,1,1,1,1)*anis_factor)) + else + f_coeff(nn,1,1,1,1)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,1))) + end if + ! + end do + else + !Save f as VX + if(RIM_id_epsm1_reference/=0) then + ! + !vX is included in the anis_factor + ! + f_coeff(1,1,1,1,1)=1.0/(4*pi*(1.0+vX_nn(1,1,1,1,1))) + call msg('r','an f_lim =',real(f_coeff(1,1,1,1,1),SP)) + else + f_coeff(1,1,1,1,1)=vX_nn(1,1,1,1,1)/(4*pi*(1.0+vX_nn(1,1,1,1,1))) + call msg('r','f_lim =',real(f_coeff(1,1,1,1,1),SP)) + end if + ! + do nn=2,4 + idm = nn-1 + ! + if(RIM_id_epsm1_reference/=0 .and. iomega==1) then + anis_factor=dot_product(em1_anis,q_grid_b_cc(:,idm)**2) + f_coeff(nn,1,1,1,1)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (f_coeff(1,1,1,1,1)*anis_factor)) + call msg('r','an f_nn =',real(f_coeff(nn,1,1,1,1),SP)) + else + f_coeff(nn,1,1,1,1)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,1))) + call msg('r','f_nn =',real(f_coeff(nn,1,1,1,1),SP)) + end if + ! + enddo + end if + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("metal") + ! + call msg("r","Metallic q --> 0 shape of f") + ! + ! The first point in metals is computed with a first order Taylor expansion + ! Coefficients in rlu (linear) + ! + ! We make average to compute the constant value, that is linked to density at Fermi. + ! Average along the three directions if 3D, only along the planar directions if 2D. + ! + if (cut_is_slab) then + f_coeff(1,1,1,1,1) = ((2._DP*f_func(idir(2)+1)-f_func(idir(2)+4))+ & +& (2._DP*f_func(idir(3)+1)-f_func(idir(3)+4)))/2._DP + call msg('r','2D Thomas-Fermi const =',-real(pi*f_coeff(1,igr,igc,iq1,iomega)*alat(idir(1)),SP)) + else + f_coeff(1,1,1,1,1) = ((2._DP*f_func(2)-f_func(5))+(2._DP*f_func(3)-f_func(6)) + & +& (2._DP*f_func(4)-f_func(7)))/3._DP + call msg('r','3D Thomas-Fermi const =',-real(4*pi**2*f_coeff(1,igr,igc,iq1,iomega),SP)) + end if + !f_coeff(1,1,1,1,1)=0._DP + ! + ! + ! quadratic dispersion to the neighbours + f_coeff(2,igr,igc,iq1,iomega) = (4._DP*f_func(2)-f_func(5)-3._DP*f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(1,1)) ! f1 + f_coeff(3,igr,igc,iq1,iomega) = (4._DP*f_func(3)-f_func(6)-3._DP*f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(2,2)) ! f2 + f_coeff(4,igr,igc,iq1,iomega) = (4._DP*f_func(4)-f_func(7)-3._DP*f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(3,3)) ! f3 + f_coeff(5,igr,igc,iq1,iomega) = (f_func(5)-2._DP*f_func(2)+f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(1,1)**2) ! f11 + f_coeff(6,igr,igc,iq1,iomega) = (f_func(6)-2._DP*f_func(3)+f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(2,2)**2) ! f22 + f_coeff(7,igr,igc,iq1,iomega) = (f_func(7)-2._DP*f_func(4)+f_coeff(1,1,1,1,1))/(2.0*q_grid_b_rlu(3,3)**2) ! f33 + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("Dirac") + ! + ! Dirac cone treatment, a Fermi velocity has to be provided + call msg('r','Graphene q--> 0 term included') + call parser('v_fermi',Dirac_v_f) + Dirac_v_f_au = Dirac_v_f/2.187691_SP + X0_Dirac = -1.0/(4._SP*Dirac_v_f_au*alat(idir(1))) + ! + call msg("r","Dirac cone q --> 0 shape of f") + call msg("r","Fermi velocity vf ",Dirac_v_f_au) + ! + !f_coeff(1,1,1,1,1) = (vX_nn(1,1,1,1,1)+VX_Dirac)/((1._DP+(vX_nn(1,1,1,1,1)+VX_Dirac)) & +!& *2._DP*pi*alat(idir(1))) + f_coeff(1,1,1,1,1) = X0_Dirac*2._DP*pi*alat(idir(1)) + ! + call msg("r","Dirac f head",real(vX_nn(1,1,1,1,1),SP)) + ! + do nn=2,3 + idm = nn-1 + ! + f_coeff(nn,1,1,1,1)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (v_norm(q_grid_b_cc(:,idm))*f_coeff(1,1,1,1,1))) + ! + enddo + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case default + ! + call msg("r","No specific q --> 0 shape of f") + ! + end select + !Do not convert static head to iku + cycle + ! + end if + ! + ! Finite frequency, as semiconductors + ! The finite frequency term should be united to the term for semiconductors, but there is to see the role of anistropy + ! + if(igr==1.and.igc==1.and.iomega.ne.1) then + select case (trim(rimw_type)) + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("semiconductor","Dirac") + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + if (cut_is_slab) then + !Save f as VXV + if(RIM_id_epsm1_reference/=0) then + ! + !vX is included in the anis_factor + ! + f_coeff(1,1,1,1,iomega)=(4._DP*pi)/(bare_qpg(1,1)*bare_qpg(1,1)*(2._DP*pi*alat(idir(1)))**2) + else + f_coeff(1,1,1,1,iomega)=vX_nn(1,1,1,1,iomega)*(4._DP*pi) / & +& (bare_qpg(1,1)*bare_qpg(1,1)*(2._DP*pi*alat(idir(1)))**2) + end if + ! + do nn=2,4 + ! + idm=nn-1 + ! + if(RIM_id_epsm1_reference/=0 .and. iomega==1) then + anis_factor=dot_product(em1_anis,q_grid_b_cc(:,idm)**2) + f_coeff(nn,1,1,1,iomega)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (f_coeff(1,1,1,1,iomega)*anis_factor)) + else + f_coeff(nn,1,1,1,iomega)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,iomega))) + end if + ! + enddo + else + !Save f as VX + f_coeff(1,1,1,1,iomega)=vX_nn(1,1,1,1,iomega)/(4*pi*(1.0+vX_nn(1,1,1,1,iomega))) + ! + do nn=2,4 + ! + idm=nn-1 + ! + if(RIM_id_epsm1_reference/=0 .and. iomega==1) then + anis_factor=dot_product(em1_anis,q_grid_b_cc(:,idm)**2) + f_coeff(nn,1,1,1,iomega)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (f_coeff(1,1,1,1,iomega)*anis_factor)) + else + f_coeff(nn,1,1,1,iomega)=1._DP/q_grid_b_rlu(idm,idm)*log(f_func(nn)/ & +& (v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,iomega))) + end if + ! + enddo + end if + cycle + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("metal") + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! + ! As in the static case, the first point in metals is missing + ! We know f behaves quadratically for q->0 + ! We fix the q->0 value from third order expansion with nnn + ! + if (cut_is_slab) then + f_coeff(1,1,1,1,iomega)= (8._DP*f_func(idir(2)+1)-f_func(idir(2)+4))/(8._DP* & +& v_norm(q_grid_b_cc(:,idir(2)))**2) & +& + (8._DP*f_func(idir(3)+1)-f_func(idir(3)+4))/(8._DP* & +& v_norm(q_grid_b_cc(:,idir(3)))**2) + else + f_coeff(1,1,1,1,iomega)= (8._DP*f_func(2)-f_func(5))/(12._DP*v_norm(q_grid_b_cc(:,1))**2) & +& + (8._DP*f_func(3)-f_func(6))/(12._DP*v_norm(q_grid_b_cc(:,2))**2) & +& + (8._DP*f_func(4)-f_func(7))/(12._DP*v_norm(q_grid_b_cc(:,3))**2) + end if + + ! Alternative expansion through third order, easier for inclusion of anistropy but more difficult to use + ! f_coeff(5,1,1,1,iomega)= (8._DP*f_func(2)-f_func(5))/(4._DP*v_norm(q_grid_b_rlu(:,1))**2) + !f_coeff(6,1,1,1,iomega)= (8._DP*f_func(3)-f_func(6))/(4._DP*v_norm(q_grid_b_rlu(:,2))**2) + !f_coeff(7,1,1,1,iomega)= (8._DP*f_func(4)-f_func(7))/(4._DP*v_norm(q_grid_b_rlu(:,3))**2) + !f_coeff(8,1,1,1,iomega)= (f_func(5)-4._DP*f_func(2))/(4._DP*v_norm(q_grid_b_rlu(:,1))**3) + !f_coeff(9,1,1,1,iomega)= (f_func(6)-4._DP*f_func(3))/(4._DP*v_norm(q_grid_b_rlu(:,2))**3) + !f_coeff(10,1,1,1,iomega)= (f_func(7)-4._DP*f_func(4))/(4._DP*v_norm(q_grid_b_rlu(:,3))**3) + do nn=2,4 + ! + idm=nn-1 + f_coeff(nn,1,1,1,iomega)= 1.0/q_grid_b_rlu(idm,idm)*log((2*f_func(nn)-sqrt(4*f_func(nn)**2 & +& -v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,iomega)* f_func(nn+3))) & +& /(2._DP*v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1,1,iomega))) + f_coeff(nn+3,1,1,1,iomega)= (4*f_func(nn)**2-f_func(nn+3)*v_norm(q_grid_b_cc(:,idm))**2*f_coeff(1,1,1, & +& 1,iomega)+2*f_func(nn)*sqrt(4*f_func(nn)**2-v_norm(q_grid_b_cc(:,idm))**2 & +& *f_coeff(1,1,1,1,iomega)*f_func(nn+3)))/ & +& (f_func(nn+3)*v_norm(q_grid_b_cc(:,idm))**2*q_grid_b_rlu(idm,idm)) + end do + ! + ! keep the coefficients in rlu + cycle + end select + ! + end if + end if + ! + ! Wings + if ( ( (igr==1) .and. (igc .ne. 1) ) .or. ( (igc==1) .and. (igr .ne. 1) ) ) then + select case (trim(rimw_type)) + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! case ("semiconductor") + ! + ! Wings for semiconductor (now removed from the calculation) + ! + !======= + ! dp_dummy(1) = f_func(1) !f0 + ! dp_dummy(2) = (-f_func(5)+4._DP*f_func(3)-3._DP*f_func(1))/(2._DP*q_grid_b_rlu(1,1)) !f1 + ! dp_dummy(3) = (-f_func(4)+4._DP*f_func(2)-3._DP*f_func(1))/(2._DP*q_grid_b_rlu(2,2)) !f2 + ! dp_dummy(4) = (f_func(5)-2._DP*f_func(3)+f_func(1))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 + ! dp_dummy(5) = (f_func(4)-2._DP*f_func(2)+f_func(1))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + case ("metal") + ! + dp_dummy(:) = 0.0 + ! + if (cut_is_slab) then + if (iomega==1) then + ! + ! Static wings in 2D metals computed with Taylor expansion, vanish linearly for q --> 0 + ! + if ( (igr==1 .and. (g_vec(igc,idir(2))==g_vec(igc,idir(3))) .and. (g_vec(igc,idir(2))==0)) .or. & +& (igc==1 .and. (g_vec(igr,idir(2))==g_vec(igr,idir(3))) .and. (g_vec(igr,idir(2))==0)) ) then + ! Out of plane wings + ! Coefficients in rlu + f_coeff(1,igr,igc,1,iomega)= 0 + dp_dummy(2)=(4.0*f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)) + dp_dummy(3)=(4.0*f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)) + dp_dummy(4)=(4.0*f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)) + dp_dummy(5)=(f_func(5)-2.0*f_func(2))/(2._DP*q_grid_b_rlu(1,1)**2) + dp_dummy(6)=(f_func(6)-2.0*f_func(3))/(2._DP*q_grid_b_rlu(2,2)**2) + dp_dummy(7)=(f_func(7)-2.0*f_func(4))/(2._DP*q_grid_b_rlu(3,3)**2) + ! linear at leading order, there is a small anistropy that is accounted averaging + + f_coeff(idir(2)+1,igr,igc,1,iomega) = (dp_dummy(idir(2)+1)+dp_dummy(idir(3)+1))/2._DP !f idir2 + f_coeff(idir(3)+1,igr,igc,1,iomega) = (dp_dummy(idir(2)+1)+dp_dummy(idir(3)+1))/2._DP !f idir3 + f_coeff(idir(1)+1,igr,igc,1,iomega) = 0 + f_coeff(idir(2)+4,igr,igc,1,iomega) = (dp_dummy(idir(2)+4)+dp_dummy(idir(3)+4))/2._DP !f idir22 + f_coeff(idir(3)+4,igr,igc,1,iomega) = (dp_dummy(idir(2)+4)+dp_dummy(idir(3)+4))/2._DP !f idir33 + f_coeff(idir(1)+4,igr,igc,1,iomega) = 0 + else + ! On plane wings + !Coefficients in rlu + f_coeff(1,igr,igc,1,iomega)= 0 + f_coeff(2,igr,igc,1,iomega) = (f_func(2)+f_func(5))/(2._DP*q_grid_b_rlu(1,1)) !f1 + f_coeff(3,igr,igc,1,iomega) = (f_func(3)+f_func(6))/(2._DP*q_grid_b_rlu(2,2)) !f2 + f_coeff(4,igr,igc,1,iomega) = (f_func(4)+f_func(7))/(2._DP*q_grid_b_rlu(3,3)) !f3 + f_coeff(5:7,igr,igc,1,iomega)= 0 + f_coeff(8,igr,igc,1,iomega) = (f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)**2) !f11a + f_coeff(9,igr,igc,1,iomega) = (f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)**2) !f22a + f_coeff(10,igr,igc,1,iomega) = (f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)**2) !f33a + ! + end if + ! keep the coefficients in rlu + cycle + else + ! + ! Wings f function, XV and W are anisotropic at finite frequency + ! Naive anisotropic expansion, split left and right neighbours + ! + ! We need the second order neighbours along the left and right direction + do nn=2,3 + iq2 = idx_q_2(igr,1,idir(nn)+1) + ig1 = idx_G_2(igr,1,idir(nn)+1) + ig2 = idx_G_2(igc,1,idir(nn)+1) + f_func(idir(nn)+7) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nnn(idir(nn)+1,1,igr,igc,iomega)/(vX_nnn(idir(nn)+1,1,igr,igc,iomega)+1._DP) + iq2 = idx_q_2(igr,1,idir(nn)+4) + ig1 = idx_G_2(igr,1,idir(nn)+4) + ig2 = idx_G_2(igc,1,idir(nn)+4) + f_func(idir(nn)+10) = real(bare_qpg(iq2,ig1)*bare_qpg(iq2,ig2),kind=DP)/(4._DP*pi)*& +& vX_nnn(idir(nn)+4,1,igr,igc,iomega)/(vX_nnn(idir(nn)+4,1,igr,igc,iomega)+1._DP) + end do + ! + if ( (igr==1 .and. (g_vec(igc,idir(2))==g_vec(igc,idir(3))) .and. (g_vec(igc,idir(2))==0)) .or. & +& (igc==1 .and. (g_vec(igr,idir(2))==g_vec(igr,idir(3))) .and. (g_vec(igr,idir(2))==0)) ) then + ! Coefficients kept in rlu + ! Out of plane plane wings quadratic + cubic (1-3 positive 5-7 negative quadratic direction) + ! (8-10 positive 11-13 negative cubic direction) + f_coeff(1:13,igr,igc,1,iomega) = 0._DP + f_coeff(idir(2)+1,igr,igc,1,iomega) = (8*f_func(idir(2)+1)-f_func(idir(2)+7)) / & +& (4*q_grid_b_rlu(idir(2),idir(2))**2) ! f idir22 r + f_coeff(idir(3)+1,igr,igc,1,iomega) = (8*f_func(idir(3)+1)-f_func(idir(3)+7)) / & +& (4*q_grid_b_rlu(idir(3),idir(3))**2) ! f idir33 r + f_coeff(idir(2)+7,igr,igc,1,iomega) = (f_func(idir(2)+7)-4*f_func(idir(2)+1)) / & +& (4*q_grid_b_rlu(idir(2),idir(2))**3) ! f idir222 r + f_coeff(idir(3)+7,igr,igc,1,iomega) = (f_func(idir(3)+7)-4*f_func(idir(3)+1)) / & + (4*q_grid_b_rlu(idir(3),idir(3))**3) ! f idir333 r + ! + f_coeff(idir(2)+4,igr,igc,1,iomega) = (8*f_func(idir(2)+4)-f_func(idir(2)+10)) / & +& (4*q_grid_b_rlu(idir(2),idir(2))**2) ! f idir22 l + f_coeff(idir(3)+4,igr,igc,1,iomega) = (8*f_func(idir(3)+4)-f_func(idir(3)+10)) / & +& (4*q_grid_b_rlu(idir(3),idir(3))**2) ! f idir33 l + f_coeff(idir(2)+10,igr,igc,1,iomega) = -(f_func(idir(2)+10)-4*f_func(idir(2)+4)) / & +& (4*q_grid_b_rlu(idir(2),idir(2))**3) ! f idir222 l + f_coeff(idir(3)+10,igr,igc,1,iomega) = -(f_func(idir(3)+10)-4*f_func(idir(3)+4)) / & +& (4*q_grid_b_rlu(idir(3),idir(3))**3) ! f idir333 l + ! + else + ! Coefficients kept in rlu + ! In plane wings linear + quadratic (1-3 positive 5-7 negative linear direction) + ! (7-9 positive 10-12 negative quadratic direction) + f_coeff(1:13,igr,igc,1,iomega) = 0._DP + f_coeff(idir(2)+1,igr,igc,1,iomega) = (4*f_func(idir(2)+1)-f_func(idir(2)+7)) / & +& (2*q_grid_b_rlu(idir(2),idir(2))) ! f idir2 r + f_coeff(idir(3)+1,igr,igc,1,iomega) = (4*f_func(idir(3)+1)-f_func(idir(3)+7)) / & +& (2*q_grid_b_rlu(idir(3),idir(3))) ! f idir3 r + f_coeff(idir(2)+7,igr,igc,1,iomega) = (f_func(idir(2)+1)-2*f_func(idir(2)+7)) / & +& (2*q_grid_b_rlu(idir(2),idir(2))**2) ! f idir22 r + f_coeff(idir(3)+7,igr,igc,1,iomega) = (f_func(idir(3)+1)-2*f_func(idir(3)+7)) / & +& (2*q_grid_b_rlu(idir(3),idir(3))**2) ! f idir33 r + ! + f_coeff(idir(2)+4,igr,igc,1,iomega) = -(4*f_func(idir(2)+4)-f_func(idir(2)+10)) / & +& (2*q_grid_b_rlu(idir(2),idir(2))) ! f idir2 l + f_coeff(idir(3)+4,igr,igc,1,iomega) = -(4*f_func(idir(3)+4)-f_func(idir(3)+10)) / & +& (2*q_grid_b_rlu(idir(3),idir(3))) ! f idir3 l + f_coeff(idir(2)+10,igr,igc,1,iomega) = (f_func(idir(2)+10)-2*f_func(idir(2)+4)) / & +& (2*q_grid_b_rlu(idir(2),idir(2))**2) ! f idir22 l + f_coeff(idir(3)+10,igr,igc,1,iomega) = (f_func(idir(3)+10)-2*f_func(idir(3)+4)) / & +& (2*q_grid_b_rlu(idir(3),idir(3))**2) ! f idir33 l + end if + ! keep the coefficients in rlu + cycle + end if + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + else + ! + ! Static wings in 3D metals computed with Taylor expansion, vanish quadratically for q --> 0 + ! Anistropy not relevant at leading order, quadratic behaviour originates from potential + ! + !======= + ! + !Coefficients in rlu + if (iomega==1) then + dp_dummy(5) = (f_func(5)+f_func(2))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 + dp_dummy(6) = (f_func(6)+f_func(3))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + dp_dummy(7) = (f_func(7)+f_func(4))/(2._DP*q_grid_b_rlu(3,3)**2) !f33 + dp_dummy(8) = (f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)**3) !f111 + dp_dummy(9) = (f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)**3) !f222 + dp_dummy(10) = (f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)**3) !f333 + else + dp_dummy(2) = (f_func(2)-f_func(5))/(2._DP*q_grid_b_rlu(1,1)) !f1 + dp_dummy(3) = (f_func(3)-f_func(6))/(2._DP*q_grid_b_rlu(2,2)) !f2 + dp_dummy(4) = (f_func(4)-f_func(7))/(2._DP*q_grid_b_rlu(3,3)) !f3 + dp_dummy(5) = (f_func(2)+f_func(5))/(2._DP*q_grid_b_rlu(1,1)**2) !f11 + dp_dummy(6) = (f_func(3)+f_func(6))/(2._DP*q_grid_b_rlu(2,2)**2) !f22 + dp_dummy(7) = (f_func(4)+f_func(7))/(2._DP*q_grid_b_rlu(3,3)**2) !f33 + + end if + ! + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + end if + ! + end select + ! + end if + end if ! + ! set to zero coefficients in the non periodic direction + if (cut_is_slab) then + do nn=0,2 + dp_dummy(1+3*nn+idir(1))=0 + end do + end if ! - !Coefficients in iku + ! Coefficients in iku + ! + f_coeff(1,igr,igc,iq1,iomega) = dp_dummy(1) !f (iku) + f_coeff(2,igr,igc,iq1,iomega) = dp_dummy(2)*a(1,1)/alat(1)+dp_dummy(3)*a(2,1)/alat(1)+dp_dummy(4)*a(3,1)/alat(1) !fx (iku) + f_coeff(3,igr,igc,iq1,iomega) = dp_dummy(3)*a(2,2)/alat(2)+dp_dummy(2)*a(1,2)/alat(2)+dp_dummy(4)*a(3,2)/alat(2) !fy (iku) + f_coeff(4,igr,igc,iq1,iomega) = dp_dummy(4)*a(3,3)/alat(3)+dp_dummy(2)*a(1,3)/alat(3)+dp_dummy(3)*a(2,3)/alat(3) !fz (iku) + !fxx(iku) + f_coeff(5,igr,igc,iq1,iomega) = dp_dummy(5)*(a(1,1)/alat(1))**2+dp_dummy(6)*(a(2,1)/alat(1))**2 + dp_dummy(7)* & +& (a(3,1)/alat(1))**2+2*dp_dummy(11)*(a(1,1)/alat(1))*(a(2,1)/alat(1)) + 2*dp_dummy(12) & +& *(a(1,1)/alat(1))*(a(3,1)/alat(1))+ 2*dp_dummy(13)*(a(2,1)/alat(1))*(a(3,1)/alat(1)) + !fxy(iku) + f_coeff(6,igr,igc,iq1,iomega) = dp_dummy(5)*a(1,1)/alat(1)*a(1,2)/alat(2)+dp_dummy(6)*a(2,2)/alat(2)*a(2,1)/alat(1) & +& +dp_dummy(7)*a(3,1)/alat(1)*a(3,2)/alat(2)+dp_dummy(11)*(a(1,1)/alat(1)*a(2,2)/alat(2) & +& +a(1,2)/alat(2)*a(2,1)/alat(1)) + dp_dummy(12)*(a(1,1)/alat(1)*a(3,2)/alat(2) & +& +a(1,2)/alat(2)*a(3,1)/alat(1)) + dp_dummy(13)*(a(3,1)/alat(1)*a(2,2)/alat(2) & +& +a(3,2)/alat(2)*a(2,1)/alat(1)) + !fxz(iku) + f_coeff(7,igr,igc,iq1,iomega) = dp_dummy(5)*a(1,1)/alat(1)*a(1,3)/alat(3)+dp_dummy(7)*a(3,3)/alat(3) *a(3,1)/alat(1) & +& +dp_dummy(6)*a(2,1)/alat(1)*a(2,3)/alat(3)+dp_dummy(11)*(a(1,1)/alat(1)*a(2,3)/alat(3) & +& +a(1,3)/alat(3)*a(2,1)/alat(1)) + dp_dummy(12)*(a(1,1)/alat(1)*a(3,3)/alat(3) & +& +a(1,3)/alat(3)*a(3,1)/alat(1)) + dp_dummy(13)*(a(2,1)/alat(1)*a(3,3)/alat(3) & +& +a(2,3)/alat(3)*a(3,1)/alat(1)) + !fyy(iku) + f_coeff(8,igr,igc,iq1,iomega) = dp_dummy(6)*(a(2,2)/alat(2))**2+dp_dummy(5)*(a(1,2)/alat(2))**2 + dp_dummy(7)* & +& (a(3,2)/alat(2))**2+2*dp_dummy(11)*(a(1,2)/alat(2))*(a(2,2)/alat(2)) + 2*dp_dummy(12) & +& *(a(1,2)/alat(2))*(a(3,2)/alat(2))+ 2*dp_dummy(13)*(a(2,2)/alat(2))*(a(3,2)/alat(2)) + !fyz(iku) + f_coeff(9,igr,igc,iq1,iomega) = dp_dummy(6)*a(2,2)/alat(2)*a(2,3)/alat(3)+dp_dummy(7)*a(3,3)/alat(3) *a(3,2)/alat(2) & +& +dp_dummy(5)*a(1,2)/alat(2)*a(1,3)/alat(3)+dp_dummy(11)*(a(1,3)/alat(3)*a(2,2)/alat(2) & +& +a(1,2)/alat(2)*a(2,3)/alat(3)) + dp_dummy(12)*(a(1,2)/alat(2)*a(3,3)/alat(3) & +& +a(1,3)/alat(3)*a(3,2)/alat(2)) + dp_dummy(13)*(a(2,2)/alat(2)*a(3,3)/alat(3) & +& +a(2,3)/alat(3)*a(3,2)/alat(2)) + !fzz(iku) + f_coeff(10,igr,igc,iq1,iomega)= dp_dummy(7)*(a(3,3)/alat(3))**2+dp_dummy(5)*(a(1,3)/alat(3))**2 + dp_dummy(6)* & +& (a(2,3)/alat(3))**2+2*dp_dummy(11)*(a(1,3)/alat(3))*(a(2,3)/alat(3)) + 2*dp_dummy(12) & +& *(a(1,3)/alat(3))*(a(3,3)/alat(3))+ 2*dp_dummy(13)*(a(2,3)/alat(3))*(a(3,3)/alat(3)) + !third order terms + f_coeff(11:20,igr,igc,iq1,iomega)=0 ! - f_coeff(1,igr,igc,iq1,iomega) = dp_dummy(1) !f (iku) - f_coeff(2,igr,igc,iq1,iomega) = dp_dummy(2)+dp_dummy(3)*a(2,1)/alat(1) !fx (iku) - f_coeff(3,igr,igc,iq1,iomega) = dp_dummy(3)+dp_dummy(2)*a(1,2)/alat(2) !fy (iku) - f_coeff(4,igr,igc,iq1,iomega) = dp_dummy(4)+dp_dummy(5)*(a(2,1)/alat(1))**2 !fxx (iku) - f_coeff(5,igr,igc,iq1,iomega) = dp_dummy(4)*a(1,2)/alat(2)+dp_dummy(5)*a(2,1)/alat(1) !fxy (iku) - f_coeff(6,igr,igc,iq1,iomega) = dp_dummy(5)+dp_dummy(4)*(a(1,2)/alat(2))**2 !fyy (iku) + do i1=1,3 !i + do i2=1,3 !j + do nn=8,10 + f_coeff(7+3*i1+i2,igr,igc,iq1,iomega) = f_coeff(7+3*i1+i2,igr,igc,iq1,iomega) + & +& dp_dummy(nn)*a(nn-7,i1)/alat(i1)*a(nn-7,i1)/alat(i1)*a(nn-7,i2)/alat(i2) !fiij(iku) + end do + end do + end do + do nn=8,10 + f_coeff(20,igr,igc,iq1,iomega) = f_coeff(20,igr,igc,iq1,iomega) + dp_dummy(nn)*a(nn-7,1)/alat(1)* & +& a(nn-7,2)/alat(2)*a(nn-7,3)/alat(3) !fxyz(iku) + end do + ! +! old version , in general not true f_coeff(:,igc,igr,iq1,iomega) = f_coeff(:,igr,igc,iq1,iomega) ! - endif - ! - f_coeff(:,igc,igr,iq1,iomega) = f_coeff(:,igr,igc,iq1,iomega) - ! +! end if + enddo enddo - enddo - enddo + enddo + !call msg('r','f111',real(f_coeff(11,1,1,2,iomega)),SP) + !call msg('r','f112',real(f_coeff(12,1,1,2,iomega)),SP) + !call msg('r','f221',real(f_coeff(14,1,1,2,iomega)),SP) + !call msg('r','f222',real(f_coeff(15,1,1,2,iomega)),SP) + !call msg('r','f123',real(f_coeff(20,1,1,2,iomega)),SP) enddo ! call timing('RIM-W-coeff',OPR='stop') @@ -276,12 +1085,16 @@ subroutine QP_interpolate_W(X,Xw,q,mode) !Calculation of RIM-W integrals !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! + call parser('rimw_check',l_check) + if(l_check) call check_W_interpolation(q,Xw,f_coeff,q_grid_b_rlu,vX_nn,idx_q,idx_is,idx_G) + ! call rim('c',Xw) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !Debugging output !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! + ! if(allocated(RIM_qpg)) then call msg('nr','Comparison between head Wc = W - v averages: < Wc > [au] RIM-W/RIM') do iq1=1,q%nibz,2 @@ -343,18 +1156,21 @@ subroutine QP_interpolate_W(X,Xw,q,mode) ! dummy(:) = f_coeff(1,1,1,1,:) if (RIM_id_epsm1_reference/=0) then - dummy(1) = dummy(1)*0.5_SP*(em1_anis(idir(2))+em1_anis(idir(3))) + if (cut_is_slab) then + dummy(1) = dummy(1)*0.5*(em1_anis(idir(2))+em1_anis(idir(3))) + else + dummy(1) = dummy(1)/3.0*(em1_anis(1)+em1_anis(2)+em1_anis(3)) + end if end if ! else ! - ! AF: a tentative implementation with a mixing constant of 1.E-5 - ! instead of 1.E-4 as below has been proposed, but led to - ! numerical differencies in the results. To be investigated + dummy(:) = f_coeff(1,ig1,ig2,1,:) & + + 1.E-5*(f_coeff(2,ig1,ig2,1,:)+f_coeff(3,ig1,ig2,1,:)+f_coeff(4,ig1,ig2,1,:) & +& + 1.E-5*(f_coeff(5,ig1,ig2,1,:)+f_coeff(8,ig1,ig2,1,:)+f_coeff(10,ig1,ig2,1,:) & +& + 2._SP*(f_coeff(6,ig1,ig2,1,:)+f_coeff(7,ig1,ig2,1,:)+f_coeff(9,ig1,ig2,1,:)))) ! - dummy(:) = f_coeff(1,ig1,ig2,1,:)+1.E-4*(f_coeff(2,ig1,ig2,1,:)+f_coeff(3,ig1,ig2,1,:)& -& +1.E-4*(f_coeff(4,ig1,ig2,1,:)+f_coeff(6,ig1,ig2,1,:)& - +2._SP*f_coeff(5,ig1,ig2,1,:))) + !call msg("r","dummy",real(dummy)) dummy(:) = dummy(:)*4._SP*pi/(bare_qpg(1,ig1)*bare_qpg(1,ig2)) dummy(:) = dummy/(1._SP-dummy) ! @@ -378,10 +1194,6 @@ subroutine QP_interpolate_W(X,Xw,q,mode) io_err=io_RIM_W(ID,mode,Xw) endif ! -#ifdef _GPU - YAMBO_ALLOC_GPU_SOURCE(DEV_VAR(RIM_W),RIM_W) -#endif - ! YAMBO_FREE(f_coeff) ! end subroutine @@ -394,29 +1206,27 @@ subroutine find_q_nns(q,q_grid_b_iku,idx_q,idx_is,idx_G) #include ! type(bz_samp):: q - real(SP) :: q_grid_b_iku(3,2) + real(SP) :: q_grid_b_iku(3,3) ! !The following variables map the extended BZ and n.n. ! - integer :: idx_q (RIM_W_ng,q%nibz,5) !q index irreduc. BZ - integer :: idx_is(RIM_W_ng,q%nibz,5) !Symm index - integer :: idx_G (RIM_W_ng,q%nibz,5) !G index - ! - !Index ordering of nearest-neighbour maps + integer :: idx_q (RIM_W_ng,q%nibz,7) !q index irreduc. BZ + integer :: idx_is(RIM_W_ng,q%nibz,7) !Symm index + integer :: idx_G (RIM_W_ng,q%nibz,7) !G index ! - ! 2 - ! ^ - ! | - ! 5 <-- 1 --> 3 - ! | - ! v - ! 4 + ! Index ordering of nearest-neighbour maps ! - ! Work Space + ! 3 4 + ! ^ / + ! |/ + ! 5 <-- 1 --> 2 + ! /| + ! L v + ! 7 6 ! integer :: iq,ig,iq_trial,ig_trial integer :: iq_trial_ibz,is_trial,i_dummy,nn_find,nn - real(SP) :: qpG(3,5),qpG_trial(3) + real(SP) :: qpG(3,7),qpG_trial(3) ! !Cycles over qpG points do iq=1,q%nibz @@ -432,10 +1242,12 @@ subroutine find_q_nns(q,q_grid_b_iku,idx_q,idx_is,idx_G) !Calculate qpG and qpG nearest neighbour qpG(:,1) = q%pt(iq,:)+g_vec(ig,:) ! - qpG(:,2) = qpG(:,1) + q_grid_b_iku(:,2) - qpG(:,3) = qpG(:,1) + q_grid_b_iku(:,1) - qpG(:,4) = qpG(:,1) - q_grid_b_iku(:,2) + qpG(:,2) = qpG(:,1) + q_grid_b_iku(:,1) + qpG(:,3) = qpG(:,1) + q_grid_b_iku(:,2) + qpG(:,4) = qpG(:,1) + q_grid_b_iku(:,3) qpG(:,5) = qpG(:,1) - q_grid_b_iku(:,1) + qpG(:,6) = qpG(:,1) - q_grid_b_iku(:,2) + qpG(:,7) = qpG(:,1) - q_grid_b_iku(:,3) ! !Find nn indexes do i_dummy=1,ng_vec @@ -454,7 +1266,7 @@ subroutine find_q_nns(q,q_grid_b_iku,idx_q,idx_is,idx_G) is_trial = q%sstar(iq_trial,2) qpG_trial = matmul(rl_sop(:,:,is_trial),q%pt(iq_trial_ibz,:)+g_vec(ig_trial,:)) ! - do nn=2,5 + do nn=2,7 if (v_is_zero(qpG_trial-qpG(:,nn))) then idx_q (ig,iq,nn) = iq_trial_ibz idx_is(ig,iq,nn) = is_trial @@ -463,18 +1275,887 @@ subroutine find_q_nns(q,q_grid_b_iku,idx_q,idx_is,idx_G) endif enddo ! - if (nn_find == 4) exit + if (nn_find == 6) exit ! enddo ! - if (nn_find == 4) exit + if (nn_find == 6) exit ! enddo ! - if(nn_find < 4) call error('Nearest neighbours searching failed in the W interpolation') + if(nn_find < 6) call error('Nearest neighbours searching failed in the W interpolation') ! enddo ! enddo ! end subroutine + ! +subroutine find_q_nnns(q,q_grid_b_iku,idx_q,idx_is,idx_G) + use pars, ONLY:SP + use R_lattice, ONLY:bz_samp,RIM_W_ng,g_vec,ng_vec,rl_sop + use vec_operate, ONLY:v_is_zero + ! +#include + ! + type(bz_samp):: q + real(SP) :: q_grid_b_iku(3,3) + ! + !The following variables map the n.n.n + ! + integer :: idx_q (RIM_W_ng,q%nibz,13) !q index irreduc. BZ + integer :: idx_is(RIM_W_ng,q%nibz,13) !Symm index + integer :: idx_G (RIM_W_ng,q%nibz,13) !G index + ! + ! Special treatment of the q-->0 term (only needed at zero frequency) + !===================================== + ! + ! Index ordering of NEXT nearest-neighbour along the three main directions + ! + ! 6 7 10 + ! ^ / / + ! | / / + ! --> 8 9 + ! ^ / ^ / + ! |/ | / + ! 2 <-- <-- 1 --> --> 5 13 1 --> (mixed direction terms) + ! /^ ^ / ^ / + ! / | | / | 12 + ! --> 11 + ! / ^ + ! / | + ! 4 3 + ! + ! Work Space + ! + integer :: iq,ig,iq_trial,ig_trial,max_numb + integer :: iq_trial_ibz,is_trial,i_dummy,nnn_find,nnn + real(SP) :: qpG(3,13),qpG_trial(3) + ! + !Cycles over qpG points + do iq=1,q%nibz + ! + do ig=1,RIM_W_ng + ! find in general 6 second order neighbours + max_numb=6 + ! For the neighbours of iq1, 12 second order neighbours are needed + if (ig==1) max_numb=12 + ! + nnn_find = 0 + ! + idx_q (ig,iq,1) = iq + idx_is(ig,iq,1) = 1 + idx_G (ig,iq,1) = ig + ! + !Calculate qpG and qpG of next nearest neighbour + qpG(:,1) = q%pt(iq,:)+g_vec(ig,:) + ! + !Calculate qpG and qpG of next nearest neighbour + qpG(:,2) = qpG(:,1) - 2.0*q_grid_b_iku(:,1) + qpG(:,3) = qpG(:,1) - 2.0*q_grid_b_iku(:,2) + qpG(:,4) = qpG(:,1) - 2.0*q_grid_b_iku(:,3) + qpG(:,5) = qpG(:,1) + 2.0*q_grid_b_iku(:,1) + qpG(:,6) = qpG(:,1) + 2.0*q_grid_b_iku(:,2) + qpG(:,7) = qpG(:,1) + 2.0*q_grid_b_iku(:,3) + if (ig==1) then + qpG(:,8) = qpG(:,1) + q_grid_b_iku(:,1) + q_grid_b_iku(:,2) + qpG(:,9) = qpG(:,1) + q_grid_b_iku(:,1) + q_grid_b_iku(:,3) + qpG(:,10) = qpG(:,1) + q_grid_b_iku(:,2) + q_grid_b_iku(:,3) + qpG(:,11) = qpG(:,1) - q_grid_b_iku(:,1) + q_grid_b_iku(:,2) + qpG(:,12) = qpG(:,1) - q_grid_b_iku(:,1) + q_grid_b_iku(:,3) + qpG(:,13) = qpG(:,1) - q_grid_b_iku(:,2) + q_grid_b_iku(:,3) + end if + !Find indexes + do i_dummy=1,ng_vec + ! + !Try first ig2 near ig + if (i_dummy < 2*ig) then + ig_trial = (2*(i_dummy-2*(i_dummy/2))-1)*(i_dummy/2) + ig + else + ig_trial = i_dummy + endif + ! + do iq_trial=1,q%nbz + ! + !Get ibz and is indexes + iq_trial_ibz = q%sstar(iq_trial,1) + is_trial = q%sstar(iq_trial,2) + qpG_trial = matmul(rl_sop(:,:,is_trial),q%pt(iq_trial_ibz,:)+g_vec(ig_trial,:)) + ! + ! nnn along positive suffice + ! + do nnn=2,max_numb+1 + if (v_is_zero(qpG_trial-qpG(:,nnn))) then + idx_q (ig,iq,nnn) = iq_trial_ibz + idx_is(ig,iq,nnn) = is_trial + idx_G (ig,iq,nnn) = ig_trial + nnn_find = nnn_find + 1 + endif + enddo + if (nnn_find == max_numb) exit + enddo + ! + if (nnn_find == max_numb) exit + ! + enddo + ! + if(nnn_find < max_numb) call error('Next nearest neighbours searching failed in the W interpolation') + ! + enddo + ! + enddo + ! +end subroutine +! +subroutine check_W_interpolation(q,Xw,f_coeff,delta_q,vX,idx_q,idx_is,idx_G) + ! + use pars, ONLY:SP,pi,zero_dfl,schlen,DP + use com, ONLY:msg + use R_lattice, ONLY:bz_samp,RIM_W_ng,g_vec,b,k_grid_b,ng_vec,rl_sop,bare_qpg,& +& RIM_W,q0_def_norm,cut_is_slab,idir,k_grid,RIM_id_epsm1_reference,RIM_epsm1,rimw_type + use vec_operate, ONLY:c2a,v_is_zero,v_norm,iku_v_norm + use X_m, ONLY:X_mat,X_t + use frequency, ONLY:w_samp + use IO_m, ONLY:OP_RD_CL,NONE,RD_CL,OP_RD,RD_CL_IF_END,& +& deliver_IO_error_message,OP_WR_CL,REP,VERIFY + use stderr, ONLY:STRING_split + use D_lattice, ONLY:alat,a + use matrix_operate, ONLY:m3inv,mat_transpose + ! + implicit none + type(bz_samp) :: q + type(w_samp) :: Xw + complex(DP) :: f_coeff(20,RIM_W_ng,RIM_W_ng,q%nibz,Xw%n_freqs) + real(SP) :: delta_q(3,3) + complex(DP) :: vX(7,q%nibz,RIM_W_ng,RIM_W_ng,Xw%n_freqs) + integer :: idx_q(RIM_W_ng,q%nibz,7),idx_G(RIM_W_ng,q%nibz,7),idx_is(RIM_W_ng,q%nibz,7) + ! + ! Work Space + ! + integer :: out_unit_W_fit, out_unit_W_num,npoints(3),find,iq,i_int,iomega + integer :: Np, i, ig1, ig2, iqibz, is, gmax, igr, igc, pos, i_dir,vecsign1,vecsign2, j + real(SP) :: dq, q_sampl(3),q_rlu(3),lcut,slab_vz1, slab_vplane1,slab_vz2,symminv(3,3),& +& slab_vplane2,pre_factor,r1,vslab,q_num_norm,q_out_norm + real(SP) :: g1_2D_mod, g2_2D_mod, qpGr(3), qpGc(3), idx_q1, q_out(3),q_num(3), & +& qpG_trial(3),anis_fact,em1_anis(3),q_min + complex(SP):: func,epsm1_sampl,W_sampl,f_coeff_loc(20) + character(schlen):: str_piece1,str_piece2,str_piece3 + character(1) :: dir(3) + integer, allocatable :: indexes(:,:,:) + integer :: q_nn, Gr_nn, Gc_nn, n_indx_steps(3) + integer :: iq_trial,iq_ibz_trial,is_trial,ig_trial + ! + Np = 51 + dq = 1./(Np-1) + q_min=1.E-5 + q_sampl = (/0.,0.,0./) + gmax = RIM_W_ng + lcut=alat(idir(1))/2._SP + dir = (/"1","2","3"/) + !out_unit_W_fit=(/200,300/) + !out_unit_W_num=(/400,500/) + out_unit_W_fit=200 + out_unit_W_num=400 + call k_ibz2bz(q,'i',.true.) + ! + !Anisotropy initialization + em1_anis=RIM_epsm1(:)-1._SP + ! + if (RIM_id_epsm1_reference<0.or.RIM_id_epsm1_reference>3) RIM_id_epsm1_reference=0 + if (RIM_id_epsm1_reference==0) em1_anis=0. + ! + n_indx_steps = k_grid/2 + !Allocate indexes + !if (delta_q(1,1) > delta_q(2,2)) then + ! n_indx_steps(1) = min(k_grid(1),k_grid(2))/2-1+1 + ! n_indx_steps(2) = max(k_grid(1),k_grid(2))/2-1+1 + !else + ! n_indx_steps(1) = max(k_grid(1),k_grid(2))/2-1+1 + ! n_indx_steps(2) = min(k_grid(1),k_grid(2))/2-1+1 + !end if + npoints = (n_indx_steps)*2+1 + allocate(indexes(max(npoints(1),npoints(2),npoints(3)),4,3)) + call msg('r', "npoints_1 = ", npoints(1)) + call msg('r', "npoints_2 = ", npoints(2)) + call msg('r', "npoints_3 = ", npoints(3)) + call msg('r', "n_indx_steps_1 = ", n_indx_steps(1)) + call msg('r', "n_indx_steps_2 = ", n_indx_steps(2)) + call msg('r', "n_indx_steps_3 = ", n_indx_steps(3)) + ! + !Print to file interpolated W data + do iomega=1,Xw%n_freqs + do i_dir=1,3 + ! + do igr=1,gmax + ! + do igc=igr,gmax + ! + !Open output files + write(str_piece1, '(I2.2)') igr + write(str_piece2, '(I2.2)') igc + write(str_piece3, '(I2.2)') iomega + open (unit=out_unit_W_fit,file="W_w="//trim(str_piece3)//"_"& +&//dir(i_dir)//"_fit_g"//trim(str_piece1)//"_g"//trim(str_piece2)//".dat",& +& action="write",status="replace") + open (unit=out_unit_W_num,file="W_w="//trim(str_piece3)//"_"& +&//dir(i_dir)//"_num_g"//trim(str_piece1)//"_g"//trim(str_piece2)//".dat",& +& action="write",status="replace") + ! + ! + !################################### + !# Get the indexes of the q points # + !################################### + ! + do i=1,npoints(i_dir) + ! + !First: get the (q+Gr,q+Gc) point + ! + q_num = (/0.,0.,0./) + q_num = delta_q(i_dir,:)*(i-1-n_indx_steps(i_dir)) + call c2a(b_in=b,v_in=q_num,mode='ka2i') + ! + qpGr = q_num + g_vec(igr,:) + qpGc = q_num + g_vec(igc,:) + ! + vecsign1 = sign(1,i-n_indx_steps(i_dir)-1) + ! + !Get q output in cc + ! + call c2a(b_in=b,v_in=q_num,mode='ki2c') + q_num_norm = v_norm(q_num)*vecsign1 + ! + ! Second: reduce (q+Gr,q+Gc) to the iBZ and store indexes + ! + find = 0 + ! + !Find nn indexes + ! + do ig_trial=1,ng_vec + ! + do iq_trial=1,q%nbz + ! + iq_ibz_trial = q%sstar(iq_trial,1) + is_trial = q%sstar(iq_trial,2) + ! + qpG_trial = matmul(rl_sop(:,:,is_trial),q%pt(iq_ibz_trial,:)+g_vec(ig_trial,:)) + ! + if (v_is_zero(qpGr-qpG_trial)) then + ! + indexes(i,1,i_dir) = iq_ibz_trial + indexes(i,2,i_dir) = is_trial + indexes(i,3,i_dir) = ig_trial + ! + find = find+1 + ! + end if + ! + if (v_is_zero(qpGc-qpG_trial)) then + ! + !Select the indexes + indexes(i,4,i_dir) = ig_trial + find = find+1 + ! + end if + ! + if (find == 2) exit + ! + end do + ! + if (find == 2) exit + ! + end do + if (find /= 2) call error('I have not find n.n. for the plotting of W') + ! Rename index variables + iq = indexes(i,1,i_dir) + is = indexes(i,2,i_dir) + ig1 = indexes(i,3,i_dir) + ig2 = indexes(i,4,i_dir) + ! + if(ig1 > RIM_W_ng.or.ig2 > RIM_W_ng) then + ! call msg("rs", "skipped ig1 =",ig1) + ! call msg("rs", "skipped ig2 =",ig2) + cycle + ! call msg("rs", "not skipped ig2") + end if + !---------------------------------- + ! Third: Evaluate numerical f,vc,W + !---------------------------------- + ! + ! Evaluate vslab + ! + vslab = 4._DP*pi/bare_qpg(iq,ig1)/bare_qpg(iq,ig2) + ! + ! Evaluate func + ! + func = bare_qpg(iq,ig1)*bare_qpg(iq,ig2)/(4._DP*pi)*vX(1,iq,ig1,ig2,iomega)/(vX(1,iq,ig1,ig2,iomega)+1._DP) + ! + ! Evaluate W and eps^-1 (for metals if ig1==ig2 save W, otherwise Wcorr) + ! + epsm1_sampl = vX(1,iq,ig1,ig2,iomega) + W_sampl = vX(1,iq,ig1,ig2,iomega)*vslab + ! + if (ig1==ig2 .and. iomega==1) then + ! + select case (trim(rimw_type)) + ! + case ("metal") + ! + epsm1_sampl = 1._SP + vX(1,iq,ig1,ig2,iomega) + ! + W_sampl = (1._SP+vX(1,iq,ig1,ig2,iomega))*vslab + ! + end select + ! + end if + ! + ! For metals, regularize the numerical limit of vc and W with Thomas-Fermi + ! + if (iq==1 .and. ig1==1 .and. ig2==1 .and. iomega==1) then + ! + select case (trim(rimw_type)) + ! + case ("metal") + ! + func = bare_qpg(iq+1,ig1)*bare_qpg(iq+1,ig2)/(4._DP*pi)*vX(1,iq+1,ig1,ig2,iomega)/(vX(1,iq+1,ig1,ig2,iomega)+1._DP) + ! + if (cut_is_slab) then + ! + epsm1_sampl = q_num_norm/(q_num_norm-2._SP*pi*func*alat(idir(1))) + ! + W_sampl = 2._SP*pi/(q_num_norm-2._SP*pi*func*alat(idir(1))) + ! + else + epsm1_sampl = q_num_norm**2/(q_num_norm**2-4._SP*pi*func) + ! + W_sampl = 4._SP*pi/(q_num_norm**2-4._SP*pi*func) + ! + end if + end select + ! + else if (iq==1 .and. ig1==1 .and. ig2.ne.1 .and. iomega==1 .and. .not. cut_is_slab) then + select case (trim(rimw_type)) + ! + case ("metal") + ! For metals, regularize the numerical limit through neighbours + ! Store f as f/V00 + func = bare_qpg(iq+1,ig2)/bare_qpg(iq+1,ig1)/(4._DP*pi)*vX(1,iq+1,ig1,ig2,iomega)/(vX(1,iq+1,ig1,ig2,iomega)+1._DP) + ! + epsm1_sampl = func/bare_qpg(iq,ig2)*bare_qpg(iq,ig1) + ! + W_sampl = (4._SP*pi)**2*func/bare_qpg(iq,ig2)**2 + ! + func = func*bare_qpg(iq,ig1)**2 + ! + end select + end if + ! + !Write output + ! + write(out_unit_W_num,"(I4, E15.6, E15.6,E15.6,E15.6, E15.6, E15.6, E15.6, E15.6, E15.6, E15.6, E15.6)") & +& i, q_num(1),q_num(2),q_num(3),q_num_norm,vslab,real(func),aimag(func),real(epsm1_sampl),& +& aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + ! + ! For first and last points, skip the fit + ! + !if(i == 1 .or. i == npoints(i_dir)) cycle + ! + !------------------------------- + ! Fourth: Evaluate fitted f,vc,W + !------------------------------- + do i_int=1,Np + ! + if (i == n_indx_steps(i_dir)+1) then + vecsign2 = sign(1,i_int-Np/2-1) + else + vecsign2=1 + end if + ! + !q_sampl in rlu + q_sampl = (/0.,0.,0./) + q_sampl = (dq*(i_int-1)-0.5)*delta_q(i_dir,:) + ! + ! keep q_sampl in rlu + q_rlu=q_sampl + ! Convert q_sampl to iku + call c2a(b_in=b,v_in=q_sampl,mode='ka2i') + ! + ! Save norm in rlu + q_out = q_sampl + call c2a(b_in=b,v_in=q_out,mode='ki2c') + q_out = q_num+q_out + q_out_norm = v_norm(q_out)*vecsign1*vecsign2 + ! + ! Evaluate v_col + !================= + r1 = 1. + if (iku_v_norm(qpGr+q_sampl) < q0_def_norm) then + r1 = r1*q0_def_norm + else + r1 = r1*iku_v_norm(qpGr+q_sampl) + end if + if (iku_v_norm(qpGc+q_sampl) < q0_def_norm) then + r1 = r1*q0_def_norm + else + r1 = r1*iku_v_norm(qpGc+q_sampl) + end if + ! + vslab=4._DP*pi + ! + if(cut_is_slab) then + !r1=iku_v_norm(qpGr+q_sampl)*iku_v_norm(qpGc+q_sampl) + !kz + slab_vz1=(q_sampl(idir(1))+qpGr(idir(1)))*2.*pi/alat(idir(1)) + slab_vz2=(q_sampl(idir(1))+qpGc(idir(1)))*2.*pi/alat(idir(1)) + !kxy + slab_vplane1=sqrt(((q_sampl(idir(2))+qpGr(idir(2)))*2.*pi/alat(idir(2)))**2+& +& ((q_sampl(idir(3))+qpGr(idir(3)))*2.*pi/alat(idir(3)))**2) + ! + slab_vplane2=sqrt(((q_sampl(idir(2))+qpGc(idir(2)))*2.*pi/alat(idir(2)))**2+& +& ((q_sampl(idir(3))+qpGc(idir(3)))*2.*pi/alat(idir(3)))**2) + ! + if ((slab_vplane1 < q0_def_norm) .and. (i_dir .ne. 3)) slab_vplane1 = q0_def_norm + if ((slab_vplane2 < q0_def_norm) .and. (i_dir .ne. 3)) slab_vplane2 = q0_def_norm + ! + pre_factor = sqrt(1.-exp(-slab_vplane1*lcut)*cos(slab_vz1*lcut))*& +& sqrt(1.-exp(-slab_vplane2*lcut)*cos(slab_vz2*lcut)) + ! + vslab=vslab*pre_factor + endif + ! + ! Evaluate func + !=============== + ! + !Transform f_coeff from iBZ to BZ + ! + call m3inv(rl_sop(:,:,is),symminv) + ! + call trans_f_coeff(f_coeff(:,ig1,ig2,iq,iomega),f_coeff_loc,symminv) + ! + func = f_coeff_loc(1) & +& +q_sampl(1)*(f_coeff_loc(2)+q_sampl(1)*f_coeff_loc(5)+2_DP*q_sampl(2)*f_coeff_loc(6)) & +& +q_sampl(2)*(f_coeff_loc(3)+q_sampl(2)*f_coeff_loc(8)+2_DP*q_sampl(3)*f_coeff_loc(9)) & +& +q_sampl(3)*(f_coeff_loc(4)+q_sampl(3)*f_coeff_loc(10)+2_DP*q_sampl(1)*f_coeff_loc(7)) & +& +q_sampl(1)**2*(q_sampl(1)*f_coeff_loc(11)+3_DP*q_sampl(2)*f_coeff_loc(12)+3_DP*q_sampl(3)*f_coeff_loc(13)) & +& +q_sampl(2)**2*(q_sampl(2)*f_coeff_loc(15)+3_DP*q_sampl(3)*f_coeff_loc(16)+3_DP*q_sampl(1)*f_coeff_loc(14)) & +& +q_sampl(3)**2*(q_sampl(3)*f_coeff_loc(19)+3_DP*q_sampl(1)*f_coeff_loc(17)+3_DP*q_sampl(2)*f_coeff_loc(18)) & +& +6*q_sampl(1)*q_sampl(2)*q_sampl(3)*f_coeff_loc(20) + ! + if(iq==1) then + ! + ! Head + ! + if(ig1==1.and.ig2==1) then + ! + ! Metal case, static and finite frequency + ! + select case (trim(rimw_type)) + ! + case ("metal") + ! + if (iomega==1) then + ! + q_out = q_sampl + call c2a(b_in=b,v_in=q_out,mode='ki2c') + func = f_coeff_loc(1) + f_coeff(2,ig1,ig2,iq,iomega)*abs(q_rlu(1)) & +& + f_coeff(3,ig1,ig2,iq,iomega)*abs(q_rlu(2))+f_coeff(4,ig1,ig2,iq,iomega)*abs(q_rlu(3)) & +& + f_coeff(5,ig1,ig2,iq,iomega)*(q_rlu(1))**2 + f_coeff(6,ig1,ig2,iq,iomega)*(q_rlu(2))**2 & +& + f_coeff(7,ig1,ig2,iq,iomega)*(q_rlu(3))**2 + ! + epsm1_sampl = 1._SP/(1._SP-vslab/r1*func) + W_sampl= vslab/r1/(1._SP-vslab/r1*func) + ! + if (iku_v_norm(q_sampl) < q0_def_norm) then + if(cut_is_slab) then + ! + epsm1_sampl = q_out_norm/(q_out_norm-2._SP*pi*func*alat(idir(1))) + ! + W_sampl = 2._SP*pi*a(idir(1),idir(1))/(q_out_norm-2._SP*pi*func*alat(idir(1))) + ! + else + ! + epsm1_sampl = q_out_norm**2/(q_out_norm**2-4._SP*pi*func) + ! + W_sampl = 4._SP*pi/(q_out_norm**2-4._SP*pi*func) + ! + end if + end if + ! + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1, real(func),& +& aimag(func),real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + cycle + else + ! + func = f_coeff(1,1,1,1,iomega) + ! + if (iku_v_norm(q_sampl) < q0_def_norm) then + ! + epsm1_sampl = vslab*func/(1.0-vslab*func) + W_sampl = func/(1.0-vslab*func)*vslab**2/q0_def_norm**2 + func= f_coeff(1,1,1,1,iomega)*iku_v_norm(q_sampl)**2 + !func= func*iku_v_norm(q_sampl)**2 + ! + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1, real(func),& +& aimag(func),real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + cycle + end if + ! + func = (func+f_coeff(5,1,1,1,iomega)*abs(q_rlu(1)) + f_coeff(6,1,1,1,iomega)*abs(q_rlu(2)) + & +& f_coeff(7,1,1,1,iomega)*abs(q_rlu(3)))*exp(sign(1._DP,real(f_coeff(3,1,1,1,iomega))) & +& *sqrt(f_coeff(2,1,1,1,iomega)**2*(q_rlu(1))**2 + f_coeff(3,1,1,1,iomega)**2*(q_rlu(2))**2 + & +& f_coeff(4,1,1,1,iomega)**2*(q_rlu(3))**2 )) + ! + func = func*iku_v_norm(q_sampl)**2 + !func = func + epsm1_sampl = func*vslab/r1/(1._SP-vslab/r1*func) + W_sampl= vslab/r1*func*vslab/r1/(1._SP-vslab/r1*func) + + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1, real(func),& +& aimag(func),real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + cycle + end if + end select + ! + ! Semiconductors and Dirac cone, static case and finite frequency + ! + if (iku_v_norm(q_sampl) < 1.e-5) then + ! + epsm1_sampl = vslab*f_coeff(1,1,1,1,iomega)/(1.0-vslab*f_coeff(1,1,1,1,iomega)) + W_sampl = f_coeff(1,1,1,1,iomega)/(1.0-vslab*f_coeff(1,1,1,1,iomega))*vslab**2/q0_def_norm**2 + if (iomega==1) then + select case (trim(rimw_type)) + case ("Dirac") + epsm1_sampl = vslab/q0_def_norm*f_coeff(1,1,1,1,iomega)/(1.0-vslab/q0_def_norm*f_coeff(1,1,1,1,iomega)) + W_sampl = f_coeff(1,1,1,1,iomega)/(1.0-vslab/q0_def_norm*f_coeff(1,1,1,1,iomega))*vslab**2/q0_def_norm**3 + end select + end if + ! + if(RIM_id_epsm1_reference /= 0 .and. iomega == 1) then + W_sampl = W_sampl*em1_anis(i_dir) + end if + ! + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1,0.0,& +& 0.0,real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + cycle + end if + ! + if (cut_is_slab) then + ! + func = f_coeff(1,1,1,1,iomega)*exp(sign(1._DP,real(f_coeff(2,1,1,1,iomega),DP))*sqrt( & +& f_coeff(2,1,1,1,iomega)**2*(q_rlu(1))**2+f_coeff(3,1,1,1,iomega)**2*(q_rlu(2))**2)) + ! + else + ! + func = f_coeff(1,1,1,1,iomega)*exp(sign(1._DP,real(f_coeff(2,1,1,1,iomega),DP))*sqrt(f_coeff(2,1,1,1,iomega)**2 & +& *(q_rlu(1))**2 + f_coeff(3,1,1,1,iomega)**2*(q_rlu(2))**2 + f_coeff(4,1,1,1,iomega)**2*(q_rlu(3))**2)) + + ! + end if + ! + if(RIM_id_epsm1_reference /= 0 .and. iomega == 1) then + anis_fact = dot_product(em1_anis,(2._SP*pi*q_sampl/alat)**2)/iku_v_norm(q_sampl)**2 + func = func*anis_fact + end if + ! + if (iomega==1) then + select case (trim(rimw_type)) + ! + case ("semiconductor") + call msg("rs", "Entered in the head check mode") + func = func*iku_v_norm(q_sampl)**2 + ! + case ("Dirac") + call msg("rs", "Entered in the head check mode") + func = func*iku_v_norm(q_sampl) + ! + case default + func = f_coeff_loc(1) +q_sampl(1)*(f_coeff_loc(2)+q_sampl(1)*f_coeff_loc(5) + 2_DP*q_sampl(2) & +& *f_coeff_loc(6)) +q_sampl(2)*(f_coeff_loc(3)+q_sampl(2)*f_coeff_loc(8) + 2_DP*q_sampl(3) & +& *f_coeff_loc(9))+q_sampl(3)*(f_coeff_loc(4)+q_sampl(3)*f_coeff_loc(10) + 2_DP*q_sampl(1) & +& *f_coeff_loc(7)) + end select + ! + else + select case (trim(rimw_type)) + ! + case ("semiconductor","Dirac") + func = func*iku_v_norm(q_sampl)**2 + ! + case default + func = f_coeff_loc(1) +q_sampl(1)*(f_coeff_loc(2)+q_sampl(1)*f_coeff_loc(5) + 2_DP*q_sampl(2) & +& *f_coeff_loc(6)) +q_sampl(2)*(f_coeff_loc(3)+q_sampl(2)*f_coeff_loc(8) + 2_DP*q_sampl(3) & +& *f_coeff_loc(9))+q_sampl(3)*(f_coeff_loc(4)+q_sampl(3)*f_coeff_loc(10) + 2_DP*q_sampl(1) & +& *f_coeff_loc(7)) + end select + ! + end if + end if + ! + ! wings + ! + if(ig1==1 .and. ig2 .ne. 1) then + ! + select case (trim(rimw_type)) + ! + case ("metal") + ! + if (cut_is_slab) then + if (iomega ==1) then + func = real(f_coeff(2,1,ig2,iq,iomega))*abs(q_rlu(1)) + real(f_coeff(3,1,ig2,iq,iomega))*abs(q_rlu(2)) & +& +real(f_coeff(4,1,ig2,iq,iomega))*abs(q_rlu(3)) + aimag(f_coeff(2,1,ig2,iq,iomega))*q_rlu(1) & +& + aimag(f_coeff(3,1,ig2,iq,iomega))*q_rlu(2) + aimag(f_coeff(4,1,ig2,iq,iomega))*q_rlu(3) & +& + f_coeff(5,1,ig2,iq,iomega)*(q_rlu(1))**2 + f_coeff(6,1,ig2,iq,iomega)*(q_rlu(2))**2 & +& + f_coeff(7,1,ig2,iq,iomega)*(q_rlu(3))**2 + f_coeff(8,1,ig2,iq,iomega)*sign(q_rlu(1)**2,q_rlu(1)) & +& + f_coeff(9,1,ig2,iq,iomega)*sign(q_rlu(2)**2,q_rlu(2)) + f_coeff(10,1,ig2,iq,iomega) & +& * sign(q_rlu(3)**2,q_rlu(3)) + else + if ((g_vec(ig2,idir(2))==g_vec(ig2,idir(3))) .and. (g_vec(ig2,idir(2))==0)) then + func = ((f_coeff(2,1,ig2,iq,iomega)+f_coeff(5,1,ig2,iq,iomega)+ (f_coeff(2,1,ig2,iq,iomega) & +& -f_coeff(5,1,ig2,iq,iomega))*sign(1._SP,q_rlu(1)))*(q_rlu(1))**2 + & +& (f_coeff(3,1,ig2,iq,iomega)+f_coeff(6,1,ig2,iq,iomega)+ (f_coeff(3,1,ig2,iq,iomega) & +& -f_coeff(6,1,ig2,iq,iomega))*sign(1._SP,q_rlu(2)))*(q_rlu(2))**2 + & +& (f_coeff(4,1,ig2,iq,iomega)+f_coeff(7,1,ig2,iq,iomega)+ (f_coeff(4,1,ig2,iq,iomega) & +& -f_coeff(7,1,ig2,iq,iomega))*sign(1._SP,q_rlu(3)))*(q_rlu(3))**2 + & +& (f_coeff(8,1,ig2,iq,iomega)+f_coeff(11,1,ig2,iq,iomega)+ (f_coeff(8,1,ig2,iq,iomega) & +& -f_coeff(11,1,ig2,iq,iomega))*sign(1._SP,q_rlu(1)))*(q_rlu(1))**3 + & +& (f_coeff(9,1,ig2,iq,iomega)+f_coeff(12,1,ig2,iq,iomega)+ (f_coeff(9,1,ig2,iq,iomega) & +& -f_coeff(12,1,ig2,iq,iomega))*sign(1._SP,q_rlu(2)))*(q_rlu(2))**3 + & +& (f_coeff(10,1,ig2,iq,iomega)+f_coeff(13,1,ig2,iq,iomega)+ (f_coeff(10,1,ig2,iq,iomega) & +& -f_coeff(13,1,ig2,iq,iomega))*sign(1._SP,q_rlu(3)))*(q_rlu(3))**3)/2.0 + else + func = ((f_coeff(2,1,ig2,iq,iomega)+f_coeff(5,1,ig2,iq,iomega)+ (f_coeff(2,1,ig2,iq,iomega) & +& -f_coeff(5,1,ig2,iq,iomega))*sign(1._SP,q_rlu(1)))*(q_rlu(1)) + & +& (f_coeff(3,1,ig2,iq,iomega)+f_coeff(6,1,ig2,iq,iomega)+ (f_coeff(3,1,ig2,iq,iomega) & +& -f_coeff(6,1,ig2,iq,iomega))*sign(1._SP,q_rlu(2)))*(q_rlu(2)) + & +& (f_coeff(4,1,ig2,iq,iomega)+f_coeff(7,1,ig2,iq,iomega)+ (f_coeff(4,1,ig2,iq,iomega) & +& -f_coeff(7,1,ig2,iq,iomega))*sign(1._SP,q_rlu(3)))*(q_rlu(3)) + & +& (f_coeff(8,1,ig2,iq,iomega)+f_coeff(11,1,ig2,iq,iomega)+ (f_coeff(8,1,ig2,iq,iomega) & +& -f_coeff(11,1,ig2,iq,iomega))*sign(1._SP,q_rlu(1)))*(q_rlu(1))**2 + & +& (f_coeff(9,1,ig2,iq,iomega)+f_coeff(12,1,ig2,iq,iomega)+ (f_coeff(9,1,ig2,iq,iomega) & +& -f_coeff(12,1,ig2,iq,iomega))*sign(1._SP,q_rlu(2)))*(q_rlu(2))**2 + & +& (f_coeff(10,1,ig2,iq,iomega)+f_coeff(13,1,ig2,iq,iomega)+ (f_coeff(10,1,ig2,iq,iomega) & +& -f_coeff(13,1,ig2,iq,iomega))*sign(1._SP,q_rlu(3)))*(q_rlu(3))**2)/2.0 + end if + end if + end if + ! + if (iku_v_norm(q_sampl) < q0_def_norm) then + ! + if (cut_is_slab) then + if (iomega ==1) then + ! System 2D, f func behaves as fcoeff*q + func = (real(f_coeff(2,1,ig2,iq,iomega))*abs(q_rlu(1)) + real(f_coeff(3,1,ig2,iq,iomega))*abs(q_rlu(2)) & +& + real(f_coeff(4,1,ig2,iq,iomega))*abs(q_rlu(3)) + aimag(f_coeff(2,1,ig2,iq,iomega))*q_rlu(1) & +& + aimag(f_coeff(3,1,ig2,iq,iomega))*q_rlu(2) + aimag(f_coeff(4,1,ig2,iq,iomega))*q_rlu(3)) & +& /iku_v_norm(q_sampl) + ! + W_sampl = func/iku_v_norm(qpGc+q_sampl)**2*vslab**2/q0_def_norm + epsm1_sampl = func/iku_v_norm(qpGc+q_sampl)*vslab*sqrt(iku_v_norm(q_sampl)/q0_def_norm) + func=func*iku_v_norm(q_sampl) + else + if ( (iomega .ne. 1) .and. (g_vec(ig2,idir(2))==g_vec(ig2,idir(3))) .and. (g_vec(ig2,idir(2))==0)) then + W_sampl = 0 + epsm1_sampl = 0 + func=0 + else + func = ((f_coeff(2,1,ig2,iq,iomega)-f_coeff(5,1,ig2,iq,iomega))*abs(q_rlu(1)) + & +& +(f_coeff(3,1,ig2,iq,iomega)-f_coeff(6,1,ig2,iq,iomega))*abs(q_rlu(2)) + & +& +(f_coeff(4,1,ig2,iq,iomega)-f_coeff(7,1,ig2,iq,iomega))*abs(q_rlu(3))) & + /iku_v_norm(q_sampl)/2.0 + W_sampl = func/iku_v_norm(qpGc+q_sampl)**2*vslab**2/q0_def_norm + epsm1_sampl = func/iku_v_norm(qpGc+q_sampl)*vslab*sqrt(iku_v_norm(q_sampl)/q0_def_norm) + func=func*iku_v_norm(q_sampl) + end if + end if + ! + else + ! System 3D, the static f func behaves as fcoeff*q^2 + ! Dinamically, f func has a linear term that however cancels out due to symmetry + func = (q_sampl(1)*(q_sampl(1)*f_coeff_loc(5) + 2_SP*q_sampl(2)*f_coeff_loc(6)) & +& +q_sampl(2)*(q_sampl(2)*f_coeff_loc(8) + 2_SP*q_sampl(3)*f_coeff_loc(9)) & +& +q_sampl(3)*(q_sampl(3)*f_coeff_loc(10) + 2_SP*q_sampl(1)*f_coeff_loc(7)))/iku_v_norm(q_sampl)**2 + ! + W_sampl = func/iku_v_norm(qpGc+q_sampl)**2*vslab**2 + epsm1_sampl = func/iku_v_norm(qpGc+q_sampl)*vslab*iku_v_norm(q_sampl) + func=func*iku_v_norm(q_sampl)**2 + end if + ! + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1, real(func),& +& aimag(func),real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + cycle + end if + end select + end if + end if + ! + epsm1_sampl = func*vslab/r1/(1._SP-vslab/r1*func) + W_sampl= vslab/r1*func*vslab/r1/(1._SP-vslab/r1*func) + if (ig1==ig2 .and. iomega==1) then + select case (trim(rimw_type)) + ! + case ("metal") + ! + epsm1_sampl = 1._SP/(1._SP-vslab/r1*func) + W_sampl= vslab/r1/(1._SP-vslab/r1*func) + ! + end select + end if + ! + !if(igr==1.and.igc == 1.and.iomega==1.and.i_int==1.and.i==1) then + ! write(*,*) "f_coeff_loc" + ! write(*,*) f_coeff_loc + !end if + !if(igr==1.and.igc == 1.and.iomega==1.and.i_int==1.and.i==1) then + ! write(*,*) "f_coeff" + ! write(*,*) f_coeff(:,ig1,ig2,iq,iomega) + ! call msg("rs","Wrong point") + ! call msg("rs","iq = ", iq) + ! call msg("rs","is = ", is) + ! call msg("rs","ig1 = ", ig1) + ! call msg("rs","ig2 = ", ig2) + ! call msg("rs","qpG = ",qpGr) + ! call msg("rs","nearest neighbours") + ! call msg("rs","iq = ", idx_q(ig1,iq,2)) + ! call msg("rs","is = ", idx_is(ig1,iq,2)) + ! call msg("rs","ig = ", idx_G(ig1,iq,2)) + ! call msg("rs","Re[vX] = ", real(vX(2,iq,ig1,ig2,iomega))) + ! call msg("rs","Im[vX] = ", aimag(vX(2,iq,ig1,ig2,iomega))) + ! call msg("rs","--------------------") + ! call msg("rs","iq = ", idx_q(ig1,iq,3)) + ! call msg("rs","is = ", idx_is(ig1,iq,3)) + ! call msg("rs","ig = ", idx_G(ig1,iq,3)) + ! call msg("rs","Re[vX] = ", real(vX(3,iq,ig1,ig2,iomega))) + ! call msg("rs","Im[vX] = ", aimag(vX(3,iq,ig1,ig2,iomega))) + ! call msg("rs","--------------------") + ! call msg("rs","iq = ", idx_q(ig1,iq,4)) + ! call msg("rs","is = ", idx_is(ig1,iq,4)) + ! call msg("rs","ig = ", idx_G(ig1,iq,4)) + ! call msg("rs","Re[vX] = ", real(vX(4,iq,ig1,ig2,iomega))) + ! call msg("rs","Im[vX] = ", aimag(vX(4,iq,ig1,ig2,iomega))) + ! call msg("rs","--------------------") + ! call msg("rs","iq = ", idx_q(ig1,iq,5)) + ! call msg("rs","is = ", idx_is(ig1,iq,5)) + ! call msg("rs","ig = ", idx_G(ig1,iq,5)) + ! call msg("rs","Re[vX] = ", real(vX(5,iq,ig1,ig2,iomega))) + ! call msg("rs","Im[vX] = ", aimag(vX(5,iq,ig1,ig2,iomega))) + !end if + ! + !Convert q_sampl to iku + !call c2a(b_in=b,v_in=q_sampl,mode='ka2i') + ! + ! + ! + write(out_unit_W_fit,"(I5,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6,E15.6)") & +& i, q_out(1),q_out(2),q_out(3), q_out_norm, vslab/r1, real(func),& +& aimag(func),real(epsm1_sampl),aimag(epsm1_sampl),real(W_sampl),aimag(W_sampl) + end do + ! + end do + ! + !Close output files + ! + close(unit=out_unit_W_fit) + close(unit=out_unit_W_num) + ! + end do + end do + end do + enddo +end subroutine +! +subroutine trans_f_coeff(f_func,f_func_loc,symi) + use pars, ONLY:SP,zero_dfl,schlen,DP + ! +#include + ! + ! Routine to perform the computation of f_coeff in BZ from iBZ + ! + complex(DP) :: f_func(20) + complex(SP) :: f_func_loc(20) + real(SP) :: symi(3,3) + integer :: i1,i2,i3,i4,i5,nn + ! + f_func_loc(1) = f_func(1) + f_func_loc(2) = symi(1,1)*f_func(2)+symi(2,1)*f_func(3)+symi(3,1)*f_func(4) + f_func_loc(3) = symi(1,2)*f_func(2)+symi(2,2)*f_func(3)+symi(3,2)*f_func(4) + f_func_loc(4) = symi(1,3)*f_func(2)+symi(2,3)*f_func(3)+symi(3,3)*f_func(4) + ! + f_func_loc(5) = symi(1,1)*symi(1,1)*f_func(5)+symi(1,1)*symi(2,1)*f_func(6)*2 & + & + symi(1,1)*symi(3,1)*f_func(7)*2+symi(2,1)*symi(2,1)*f_func(8) & + & + symi(2,1)*symi(3,1)*f_func(9)*2+symi(3,1)*symi(3,1)*f_func(10) + ! + f_func_loc(6) = symi(1,1)*symi(1,2)*f_func(5)+(symi(1,1)*symi(2,2)+symi(2,1)*symi(1,2))*f_func(6) & + & + (symi(1,1)*symi(3,2)+symi(3,1)*symi(1,2))*f_func(7)+symi(2,1)*symi(2,2)*f_func(8) & + & + (symi(2,1)*symi(3,2)+symi(3,1)*symi(2,2))*f_func(9)+symi(3,1)*symi(3,2)*f_func(10) + ! + f_func_loc(7) = symi(1,1)*symi(1,3)*f_func(5)+(symi(1,1)*symi(2,3)+symi(2,1)*symi(1,3))*f_func(6) & + & + (symi(1,1)*symi(3,3)+symi(3,1)*symi(1,3))*f_func(7)+symi(2,1)*symi(2,3)*f_func(8) & + & + (symi(2,1)*symi(3,3)+symi(3,1)*symi(2,3))*f_func(9)+symi(3,1)*symi(3,3)*f_func(10) + ! + f_func_loc(8) = symi(1,2)*symi(1,2)*f_func(5)+symi(1,2)*symi(2,2)*f_func(6)*2 & + & + symi(1,2)*symi(3,2)*f_func(7)*2+symi(2,2)*symi(2,2)*f_func(8) & + & + symi(2,2)*symi(3,2)*f_func(9)*2+symi(3,2)*symi(3,2)*f_func(10) + ! + f_func_loc(9) = symi(1,2)*symi(1,3)*f_func(5)+(symi(1,2)*symi(2,3)+symi(2,2)*symi(1,3))*f_func(6) & + & + (symi(1,2)*symi(3,3)+symi(3,2)*symi(1,3))*f_func(7)+symi(2,2)*symi(2,3)*f_func(8) & + & + (symi(2,2)*symi(3,3)+symi(3,2)*symi(2,3))*f_func(9)+symi(3,2)*symi(3,3)*f_func(10) + ! + f_func_loc(10)= symi(1,3)*symi(1,3)*f_func(5)+symi(1,3)*symi(2,3)*f_func(6)*2 & + & + symi(1,3)*symi(3,3)*f_func(7)*2+symi(2,3)*symi(2,3)*f_func(8) & + & + symi(2,3)*symi(3,3)*f_func(9)*2+symi(3,3)*symi(3,3)*f_func(10) + ! + f_func_loc(11:20)=0 + do i1=1,3 !i + do i2=1,3 !j + do nn=11,19 + if (mod(nn-11,4)==0) then + i3=(nn-11)/4+1 + f_func_loc(7+3*i1+i2) = f_func_loc(7+3*i1+i2) + f_func(nn)*symi(i3,i1)*symi(i3,i1)*symi(i3,i2) + else + i4=mod(nn-11,3)+1 + i3=(nn-10-i4)/3+1 + f_func_loc(7+3*i1+i2) = f_func_loc(7+3*i1+i2) + f_func(nn)*(2.0*symi(i4,i1)*symi(i3,i1)*symi(i3,i2) & +& + symi(i4,i2)*symi(i3,i1)*symi(i3,i1)) + end if + end do + do i3=0,5 + i4=(i3-mod(i3,2))/2+2+mod(i3,2) + i5=i4+1-2*mod(i3,2) + if (i4>3) i4=i4-3 + if (i5>3) i5=i5-3 + f_func_loc(7+3*i1+i2) = f_func_loc(7+3*i1+i2) + f_func(20)*symi((i3-mod(i3,2))/2+1,i1)*symi(i4,i1)*symi(i5,i2) + end do + end do + end do + ! + do nn=11,19 + if (mod(nn-11,4)==0) then + i3=(nn-11)/4+1 + f_func_loc(20) = f_func_loc(20) + f_func(nn)*symi(i3,1)*symi(i3,2)*symi(i3,3) + else + i4=mod(nn-11,3)+1 + i3=(nn-10-i4)/3+1 + f_func_loc(20) = f_func_loc(20) + f_func(nn)*(symi(i3,1)*symi(i3,2)*symi(i4,3) + symi(i3,1)*symi(i4,2)*symi(i3,3) & +& + symi(i4,1)*symi(i3,1)*symi(i3,3)) + end if + end do + do i3=0,5 + i4=(i3-mod(i3,2))/2+2+mod(i3,2) + i5=i4+1-2*mod(i3,2) + if (i4>3) i4=i4-3 + if (i5>3) i5=i5-3 + f_func_loc(20) = f_func_loc(20) + f_func(20)*symi((i3-mod(i3,2))/2+1,1)*symi(i4,2)*symi(i5,3) + end do +end subroutine + + + + + + + + diff --git a/src/qp/QP_mpa.F b/src/qp/QP_mpa.F index beae81406d..4bbc2f968c 100644 --- a/src/qp/QP_mpa.F +++ b/src/qp/QP_mpa.F @@ -18,7 +18,7 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) use electrons, ONLY:levels,spin_occ,spin use LIVE_t, ONLY:live_timing use com, ONLY:msg - use drivers, ONLY:l_sc_run,l_RIM_W + use drivers, ONLY:l_sc_run,l_rim_w use parallel_int, ONLY:PP_wait,PP_redux_wait,PARALLEL_global_indexes,PARALLEL_WF_index,& & PARALLEL_WF_distribute,PARALLEL_index use parser_m, ONLY:parser @@ -31,11 +31,11 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) use IO_m, ONLY:manage_action,OP_RD_CL,OP_WR_CL,OP_APP_CL,REP,VERIFY,NONE,RD_CL,OP_RD,RD_CL_IF_END,& & io_RESPONSE,io_MULTIPOLE,deliver_IO_error_message use QP_m, ONLY:QP_t,QP_n_G_bands,QP_Sc_steps,QP_dSc_steps,QP_solver,& -& QP_Sc,QP_n_states,QP_G_damp,QP_table,QP_dSc_delta +& QP_Sc,QP_n_states,QP_G_damp,QP_table,QP_dSc_delta,QP_G_mprop use ALLOC, ONLY:X_ALLOC_elemental use X_m, ONLY:X_par,X_t use wave_func, ONLY:WF - use R_lattice, ONLY:qindx_S,bz_samp,RIM_W_ng,RIM_W + use R_lattice, ONLY:qindx_S,bz_samp,G_m_G,RIM_W_ng,RIM_W,RIM_qpg use D_lattice, ONLY:nsym,i_time_rev,i_space_inv,mag_syms use interfaces, ONLY:QP_state_print,WF_load,WF_free,MATRIX_transpose use timing_m, ONLY:timing @@ -59,8 +59,9 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) ! ! Work Space ! - integer :: i_qp,i_w,iqbz,iqibz,ib,ig1,ig2,iqs,i_qp_to_start,iq_to_start,is,& -& iq_mem,X_range(2),io_err,ID,IO_ACT,timing_steps + integer :: i_qp,i_w,iqbz,iqibz,ib,ig1,ig2,alloc_err,iqs,iscs_save(2,4),& +& i_qp_to_start,iq_to_start,is,iq_mem,X_range(2),io_err,ID,IO_ACT,timing_steps,& +& cond_g, cond_b, m_idx, QP_G_mprop_temp ! complex(SP), allocatable :: W_(:),dc(:) type(elemental_collision), target :: isc,iscp @@ -69,14 +70,14 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) character(schlen) :: ch,SECTION_name ! integer :: X_mpa_npoles,X_ng - logical :: X_is_TR_rotated,l_X_ALLOC_elemental,l_RIM_W_g,l_para_G_MPA + logical :: X_is_TR_rotated,l_X_ALLOC_elemental,l_rim_w_g,l_para_G_MPA real(SP) :: E_kmq,f_kmq real(DP) :: dp_dummy_r,dp_dummy_i - complex(DP) :: W_1,W_i + complex(DP) :: W_i,cmplx_phase,prefac #ifdef _GPU - complex(DP) :: dp_dummy,ctmp + complex(DP) :: dp_dummy,ctmp,poles,fract_sum #else - complex(DP), allocatable :: dp_dummy(:),ctmp(:) + complex(DP), allocatable :: dp_dummy(:),ctmp(:),poles(:),fract_sum(:) #endif ! complex(SP), pointer DEV_ATTR :: isc_rhotw_p(:) @@ -369,7 +370,7 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) endif ! - !$omp parallel default(shared),private(ig1,ig2,MPred,PPcond_rate,MP_err,cond_num,i_np,MPA_Xo,l_RIM_W_g), & + !$omp parallel default(shared),private(ig1,ig2,MPred,PPcond_rate,MP_err,cond_num,i_np,MPA_Xo,l_rim_w_g), & !$omp & reduction(+:MPred_rate,PPcond_Qrate,MP_Qerr,cond_numQ) ! YAMBO_ALLOC(MPA_Xo,(Xw%n_freqs)) @@ -378,17 +379,13 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) do ig2=X_cols(1),X_cols(2) do ig1=X_rows(1),X_rows(2) ! - l_RIM_W_g=(l_RIM_W.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) + l_rim_w_g=(l_rim_w.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) ! - if (l_RIM_W_g) then - ! + if(l_rim_w_g) then MPA_Xo(1:Xw%n_freqs)=RIM_W(1:Xw%n_freqs,iqibz,ig1,ig2)/2._SP - ! else - ! - !DALV: the multiplication by isc%gamp(ig1,ig2) is performed later - MPA_Xo(1:Xw%n_freqs)=X_par(iq_mem)%blc(ig1,ig2,1:Xw%n_freqs) - ! + ! GS the multiplication by isc%gamp(ig1,ig2) was moved later when not doing rimw + MPA_Xo(1:Xw%n_freqs)=X_par(iq_mem)%blc(ig1,ig2,1:Xw%n_freqs)*isc%gamp(ig1,ig2) endif ! MPred=.false. @@ -409,6 +406,7 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) ! enddo enddo + ! !$omp end do ! YAMBO_FREE(MPA_Xo) @@ -502,8 +500,9 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) iscp%is=(/QP_table(i_qp,2),QP_table(i_qp,3),1,spin(QP_table(i_qp,:))/) iscp%qs=isc%qs ! - ! DALV: here the grid is centered in E0 - forall (i_w=1:QP_dSc_steps) W_(i_w)=Sc_W(i_qp)%p(i_w)+cI*QP_G_damp + ! DALV: here the grid is center in E0 + ! GS: Damping of G used later directly when computing the self-energy + forall (i_w=1:QP_dSc_steps) W_(i_w)=Sc_W(i_qp)%p(i_w) ! do ib=QP_n_G_bands(1),QP_n_G_bands(2) ! @@ -528,6 +527,12 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) E_kmq=E%E(isc%os(1),isc%os(2),isc%os(4)) f_kmq=E%f(isc%os(1),isc%os(2),isc%os(4)) ! + ! order of the Lorentzian propagator, default is 1 + QP_G_mprop_temp=1 + if ((QP_G_mprop>1) .and. (abs(E_kmq) >= abs(cos(pi/(2._DP*QP_G_mprop))*QP_G_damp))) then + QP_G_mprop_temp=QP_G_mprop + end if + ! #if defined _GPU ! do i_w=1,QP_dSc_steps @@ -535,16 +540,15 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) dp_dummy_r = 0.0_DP dp_dummy_i = 0.0_DP ! - W_1=W_(1) W_i=W_(i_w) ! - !DEV_ACC_DEBUG data present(MPA_E_par_p,MPA_R_par_p,isc_rhotw_p,iscp_rhotw_p,isc_gamp_p) - !DEV_ACC parallel loop collapse(3) private(i_np,ig1,ig2,bose_PPA_E,ctmp,l_RIM_W_g) & + !DEV_ACC_DEBUG data present(MPA_E_par_p,MPA_R_par_p,isc_rhotw_p,iscp_rhotw_p) + !DEV_ACC parallel loop collapse(3) private(i_np,ig1,ig2,bose_PPA_E,ctmp,cmplx_phase,prefac,poles,fract_sum,m_idx) & !DEV_ACC reduction(+:dp_dummy_r,dp_dummy_i) !DEV_CUF kernel do(3) - !DEV_OMPGPU target map(present,alloc:MPA_E_par_p,MPA_R_par_p,isc_gamp_p,isc_rhotw_p,iscp_rhotw_p) & + !DEV_OMPGPU target map(present,alloc:MPA_E_par_p,MPA_R_par_p,isc_rhotw_p,iscp_rhotw_p) & !DEV_OMPGPU & map(tofrom:dp_dummy_r,dp_dummy_i) - !DEV_OMPGPU teams loop collapse(3) private(i_np,ig1,ig2,bose_PPA_E,ctmp,l_RIM_W_g) & + !DEV_OMPGPU teams loop collapse(3) private(i_np,ig1,ig2,bose_PPA_E,ctmp,cmplx_phase,prefac,poles,fract_sum,m_idx) & !DEV_OMPGPU & reduction(+:dp_dummy_r,dp_dummy_i) ! do i_np=1,X_mpa_npoles @@ -568,28 +572,29 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) endif endif ! - l_RIM_W_g=(l_RIM_W.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) + ! DALV: the factor isc%gamp(ig1,ig2) was included here ! - if (l_RIM_W_g) then - ! - ctmp = -4._DP/spin_occ*pi*isc_rhotw_p(ig1)*conjg(iscp_rhotw_p(ig2)) *& -& (sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *MPA_R_par_p(ig1,ig2,i_np))*& -& ( (spin_occ-f_kmq+bose_PPA_E)/(W_1-E_kmq & -& +sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np))+& -& (f_kmq+bose_PPA_E)/(conjg(W_1)-E_kmq & -& -sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) ) - ! - else - ! - ! DALV: the factor isc%gamp(ig1,ig2) is included here - ctmp = -4._DP/spin_occ*pi*isc_rhotw_p(ig1)*conjg(iscp_rhotw_p(ig2)) *isc_gamp_p(ig1,ig2) *& -& (sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *MPA_R_par_p(ig1,ig2,i_np))*& -& ( (spin_occ-f_kmq+bose_PPA_E)/(W_i-E_kmq & -& +sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np))+& -& (f_kmq+bose_PPA_E)/(conjg(W_i)-E_kmq & -& -sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) ) + ctmp=0._DP + ! + ! GS : Loop for the higher order propagator. + ! The higher order propagator is implemented as from + ! Chiarotti Marzari Ferretti, Phys. Rev. Res. 4, 013242 (2022) + ! + ! The poles with energy close to zero are always treated + ! with a first order propagator + ! + do m_idx=0,QP_G_mprop_temp-1 ! - endif + cmplx_phase = exp(cI*pi*(1._SP+2._SP*m_idx)/(2._SP*QP_G_mprop_temp)) + prefac = cI*cmplx_phase*sin(pi/2._SP*QP_G_mprop_temp) + poles = W_i+cmplx_phase*QP_G_damp-E_kmq + fract_sum = (spin_occ-f_kmq+bose_PPA_E)/ & +& (poles+sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) & +& + (f_kmq+bose_PPA_E)/ & +& (conjg(poles)-sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) + ctmp = ctmp-prefac*isc_rhotw_p(ig1)*conjg(iscp_rhotw_p(ig2))* & +& (sign(1._SP,imag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_R_par_p(ig1,ig2,i_np))*fract_sum + enddo ! dp_dummy_r=dp_dummy_r+real(ctmp,DP) dp_dummy_i=dp_dummy_i+imag(ctmp) @@ -604,7 +609,7 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) ! dp_dummy=cmplx(dp_dummy_r,dp_dummy_i,kind=DP) ! - dc(i_w) = cmplx(dp_dummy,KIND=SP) + dc(i_w) = -4._DP/spin_occ*pi*cmplx(dp_dummy,KIND=SP) ! enddo ! @@ -617,9 +622,12 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) ! dp_dummy(:) = 0.0_DP ! - !DEV_OMP parallel default(shared), private(i_np,ig1,ig2,bose_PPA_E,ctmp,l_RIM_W_g) + !DEV_OMP parallel default(shared), private(i_np,ig1,ig2,bose_PPA_E,ctmp,& + !DEV_OMP & m_idx,cmplx_phase,prefac,poles,fract_sum) ! YAMBO_ALLOC(ctmp,(QP_dSc_steps)) + YAMBO_ALLOC(poles,(QP_dSc_steps)) + YAMBO_ALLOC(fract_sum,(QP_dSc_steps)) ! !DEV_OMP do reduction(+:dp_dummy), collapse(2) ! @@ -647,29 +655,19 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) endif endif ! - l_RIM_W_g=(l_RIM_W.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) - ! - if (l_RIM_W_g) then + do m_idx=0,QP_G_mprop_temp-1 ! - ctmp(:) = ctmp(:) + & -& (sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *MPA_R_par_p(ig1,ig2,i_np)) *& -& ( (spin_occ-f_kmq+bose_PPA_E)/(W_(1)-E_kmq +sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *& -& MPA_E_par_p(ig1,ig2,i_np))+& -& (f_kmq+bose_PPA_E)/(conjg(W_(1))-E_kmq -sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *& -& MPA_E_par_p(ig1,ig2,i_np)) ) + cmplx_phase = exp(cI*pi*(1._SP+2._SP*m_idx)/(2._SP*QP_G_mprop_temp)) + prefac = cI*cmplx_phase*sin(pi/2._SP*QP_G_mprop_temp) + poles(:) = W_(:)+cmplx_phase*QP_G_damp-E_kmq + fract_sum(:)= (spin_occ-f_kmq+bose_PPA_E)/ & +& (poles(:)+sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) & +& + (f_kmq+bose_PPA_E)/ & +& (conjg(poles(:))-sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_E_par_p(ig1,ig2,i_np)) + ctmp(:) = ctmp(:) -prefac * & +& (sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np)))*MPA_R_par_p(ig1,ig2,i_np))*fract_sum(:) ! - else - ! - ! DALV: the factor isc%gamp(ig1,ig2) is included here - ctmp(:) = ctmp(:) + & -& isc_gamp_p(ig1,ig2) *& -& (sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *MPA_R_par_p(ig1,ig2,i_np)) *& -& ( (spin_occ-f_kmq+bose_PPA_E)/(W_(:)-E_kmq +sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *& -& MPA_E_par_p(ig1,ig2,i_np))+& -& (f_kmq+bose_PPA_E)/(conjg(W_(:))-E_kmq -sign(1._SP,aimag(MPA_E_par_p(ig1,ig2,i_np))) *& -& MPA_E_par_p(ig1,ig2,i_np)) ) - ! - endif + enddo ! enddo ! @@ -681,6 +679,8 @@ subroutine QP_mpa(X,Xk,E,k,q,qp,Xw,GW_iter) !DEV_OMP end do ! YAMBO_FREE(ctmp) + YAMBO_FREE(poles) + YAMBO_FREE(fract_sum) ! !DEV_OMP end parallel ! diff --git a/src/qp/QP_newton.F b/src/qp/QP_newton.F index 3507210d62..53c9cbcae8 100644 --- a/src/qp/QP_newton.F +++ b/src/qp/QP_newton.F @@ -68,7 +68,7 @@ subroutine QP_newton(X,Xen,Xk,en,k,q,qp,Xw,Dip) ! else if (l_mpa) then ! - if (l_RIM_W) call QP_interpolate_W(X,Xw,q,'MPA') + if (l_rim_w) call QP_interpolate_W(X,Xw,q,'MPA') call QP_mpa(X,Xk,en,k,q,qp,Xw,iter) ! else if (l_cohsex) then diff --git a/src/qp/QP_of.F b/src/qp/QP_of.F index 7a30ca96f6..e85490efda 100644 --- a/src/qp/QP_of.F +++ b/src/qp/QP_of.F @@ -196,10 +196,15 @@ subroutine QP_of(qp,en,QPdb_read_err,what) endif ! do i_w=1,qp%GreenF_n_steps - call OUTPUT_driver(trim(G_sc_name),TITLES=(/"Energy","IM(En)"/),& + ! + ! AF,GS: if IM(En) -> Im(En), an extra column with the damping of + ! G_1 = [ w-eps0-S(w) +/-id ]^-1 is printed and disrupt the + ! testing (currently disabled in view of this) + ! + call OUTPUT_driver(trim(G_sc_name),TITLES=(/"Energy","IM(En)"/),& R_VALUES=(/real(qp%GreenF_W(i_qp,i_w)),aimag(qp%GreenF_W(i_qp,i_w))/),UNIT="eV") call OUTPUT_driver(trim(G_sc_name),TITLES=(/"Re(G)","Im(G)"/),& -& R_VALUES=(/real(qp%GreenF(i_qp,i_w)),aimag(qp%GreenF(i_qp,i_w))/),UNIT="eVm1") +& R_VALUES=(/real(qp%GreenF(i_qp,i_w)),aimag(qp%GreenF(i_qp,i_w))/),UNIT="eVm1") if (allocated(QP_Vxc).and.allocated(QP_Vnl_xc)) then call OUTPUT_driver(trim(G_sc_name),TITLES=(/"Re(S_c)"/),& & R_VALUES=(/real(qp%S_total(i_qp,i_w)-QP_Vnl_xc(i_qp)+QP_Vxc(i_qp))/),UNIT="eV") diff --git a/src/qp/QP_ppa_cohsex.F b/src/qp/QP_ppa_cohsex.F index 38b4b47e62..99c4b9b6d2 100644 --- a/src/qp/QP_ppa_cohsex.F +++ b/src/qp/QP_ppa_cohsex.F @@ -339,17 +339,17 @@ subroutine QP_ppa_cohsex(X,Xk,E,k,q,qp,Xw,GW_iter) ! if (l_ppa) then ! - !$omp parallel do default(shared), private(ig1,ig2,l_RIM_W_g), & + !$omp parallel do default(shared), private(ig1,ig2,l_rim_w_g), & !$omp & reduction(+:PPcond_rate,TO_rate,PP_err), collapse(2) do ig2=X_par(iq_mem)%cols(1),X_par(iq_mem)%cols(2) do ig1=X_par(iq_mem)%rows(1),X_par(iq_mem)%rows(2) ! ! RIM W support ! - l_RIM_W_g=(l_RIM_W.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng.and.iqibz==1) - if (RIM_W_is_diagonal.and.l_RIM_W_g) l_RIM_W_g=(ig1==ig2) + l_rim_w_g=(l_rim_w.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng.and.iqibz==1) + if (RIM_W_is_diagonal.and.l_rim_w_g) l_rim_w_g=(ig1==ig2) ! - if (l_RIM_W_g) then + if (l_rim_w_g) then X_par(iq_mem)%blc(ig1,ig2,X_range2)=RIM_W_E(ig1,ig2) ! else @@ -699,17 +699,18 @@ subroutine QP_ppa_cohsex(X,Xk,E,k,q,qp,Xw,GW_iter) ! CUF reductions. dp_dummy_r = 0.0_DP dp_dummy_i = 0.0_DP + ! dp_dummy = 0.0_DP ! !DEV_ACC_DEBUG data present(X_blc_p,RIM_W,isc_gamp_p,conjg_iscp_rhotw,isc_rhotw_p,iscp_rhotw_p) - !DEV_ACC parallel loop collapse(2) private(PPA_E,PPA_R,l_RIM_W_g,bose_PPA_E,ctmp) & + !DEV_ACC parallel loop collapse(2) private(PPA_E,PPA_R,l_rim_w_g,bose_PPA_E,ctmp) & !DEV_ACC reduction(+:dp_dummy_r,dp_dummy_i) !DEV_CUF kernel do(2) !DEV_OMPGPU target map(present,alloc:X_blc_p,RIM_W,isc_gamp_p,conjg_iscp_rhotw,isc_rhotw_p,iscp_rhotw_p) & !DEV_OMPGPU & map(tofrom:dp_dummy_r,dp_dummy_i) - !DEV_OMPGPU teams loop collapse(2) private(PPA_E,PPA_R,l_rim_W_g,bose_PPA_E,ctmp) & + !DEV_OMPGPU teams loop collapse(2) private(PPA_E,PPA_R,l_rim_w_g,bose_PPA_E,ctmp) & !DEV_OMPGPU & reduction(+:dp_dummy_r,dp_dummy_i) - !DEV_OMP parallel do default(shared), private(ig1,ig2,PPA_E,PPA_R,ctmp,l_RIM_W_g), & + !DEV_OMP parallel do default(shared), private(ig1,ig2,PPA_E,PPA_R,ctmp,l_rim_w_g), & !DEV_OMP & reduction(+:dp_dummy), collapse(2) ! do ig2=X_cols1,X_cols2 @@ -719,10 +720,10 @@ subroutine QP_ppa_cohsex(X,Xk,E,k,q,qp,Xw,GW_iter) ! ! RIM W support ! - l_RIM_W_g=(l_RIM_W.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) - if (RIM_W_is_diagonal.and.l_RIM_W_g) l_RIM_W_g=(ig1==ig2) + l_rim_w_g=(l_rim_w.and.ig1<=RIM_W_ng.and.ig2<=RIM_W_ng) + if (RIM_W_is_diagonal.and.l_rim_w_g) l_rim_w_g=(ig1==ig2) ! - if (l_RIM_W_g) then + if (l_rim_w_g) then PPA_R=-cmplx(real(DEV_VAR(RIM_W)(1,iqibz,ig1,ig2))/2._SP, & & aimag(X_blc_p(ig1,ig2,X_range1))*real(isc_gamp_p(ig1,ig2)),kind=SP)/2._SP*PPA_E else