TABLE OF CONTENTS
- ABINIT/m_screening_driver
- m_screening_driver/calc_rpa_functional
- m_screening_driver/chi0_bksmask
- m_screening_driver/random_stopping_power
- m_screening_driver/screening
- m_screening_driver/setup_screening
ABINIT/m_screening_driver [ Modules ]
NAME
m_screening_driver
FUNCTION
Calculate screening and dielectric functions
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, GMR, VO, LR, RWG, MT, RShaltaf, AS, FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_screening_driver 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_dtset 28 use m_xmpi 29 use m_xomp 30 use m_errors 31 use m_ab7_mixing 32 use m_kxc 33 use m_nctk 34 use netcdf 35 use libxc_functionals 36 use m_hdr 37 use m_dtfil 38 use m_distribfft 39 use m_crystal 40 41 use defs_datatypes, only : pseudopotential_type, ebands_t 42 use defs_abitypes, only : MPI_type 43 use m_time, only : timab 44 use m_io_tools, only : open_file, file_exists, iomode_from_fname 45 use m_fstrings, only : int2char10, sjoin, strcat, itoa, ltoa, itoa 46 use m_energies, only : energies_type, energies_init 47 use m_numeric_tools, only : print_arr, coeffs_gausslegint, c2r 48 use m_geometry, only : normv, vdotw, mkrdim, metric 49 use m_gwdefs, only : GW_TOLQ0, GW_TOLQ, em1params_t, GW_Q0_DEFAULT 50 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 51 use m_ebands, only : ebands_update_occ, ebands_copy, ebands_get_valence_idx, ebands_get_occupied, & 52 ebands_apply_scissors, ebands_free, ebands_has_metal_scheme, ebands_ncwrite, ebands_init, & 53 gaps_t, ebands_get_gaps 54 use m_bz_mesh, only : kmesh_t, littlegroup_t, littlegroup_free, get_ng0sh, find_qmesh 55 use m_kg, only : getph 56 use m_gsphere, only : gsphere_t, setshells 57 use m_vcoul, only : vcoul_t 58 use m_qparticles, only : rdqps, rdgw, show_QP 59 use m_screening, only : make_epsm1_driver, lwl_write, chi_t, chi_free, chi_new 60 use m_io_screening, only : hscr_new, hscr_io, write_screening, hscr_t 61 use m_spectra, only : spectra_t, W_EM_LF, W_EM_NLF, W_EELF 62 use m_fftcore, only : print_ngfft 63 use m_fft_mesh, only : rotate_FFT_mesh, cigfft, get_gftt, setmesh 64 use m_fft, only : fourdp 65 use m_wfd, only : wfd_init, wfdgw_t, wfdgw_copy, test_charge 66 use m_wfk, only : wfk_read_eigenvalues 67 use m_io_kss, only : make_gvec_kss 68 use m_chi0tk, only : output_chi0sumrule 69 use m_pawang, only : pawang_type 70 use m_pawrad, only : pawrad_type 71 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 72 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 73 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 74 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 75 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,& 76 & pawrhoij_free, pawrhoij_symrhoij, pawrhoij_inquire_dim 77 use m_pawdij, only : pawdij, symdij_all 78 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 79 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 80 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 81 use m_paw_sphharm, only : setsym_ylm 82 use m_paw_onsite, only : pawnabla_init 83 use m_paw_nhat, only : nhatgrid, pawmknhat 84 use m_paw_denpot, only : pawdenpot 85 use m_paw_init, only : pawinit, paw_gencond 86 use m_paw_tools, only : chkpawovlp,pawprt 87 use m_chi0, only : cchi0, cchi0q0, chi0q0_intraband 88 use m_setvtr, only : setvtr 89 use m_mkrho, only : prtrhomxmn 90 use m_pspini, only : pspini 91 use m_paw_correlations, only : pawpuxinit 92 use m_plowannier, only : plowannier_type,init_plowannier,get_plowannier, fullbz_plowannier,destroy_plowannier 93 !#ifdef __HAVE_GREENX 94 use minimax_grids, only : gx_minimax_grid !, gx_get_error_message 95 !#endif 96 97 implicit none 98 99 private
m_screening_driver/calc_rpa_functional [ Functions ]
[ Top ] [ m_screening_driver ] [ Functions ]
NAME
calc_rpa_functional
FUNCTION
Routine used to calculate the Galitskii-Migdal and RPA approximations to the correlation energy from the irreducible polarizability.
INPUTS
iq=index of the q-point in the array Qmesh%ibz where epsilon^-1 has to be calculated Ep<em1params_t>=Structure with parameters and dimensions related to the inverse dielectric matrix. Pvc<vcoul_t>=Structure gathering data on the Coulombian interaction Qmesh<kmesh_t>=Data type with information on the q-sampling Dtfil<Datafiles_type)>=variables related to files spaceComm=MPI communicator.
OUTPUT
SOURCE
2600 subroutine calc_rpa_functional(gwrpacorr,gwgmcorr,iqcalc,iq,Ep,Pvc,Qmesh,Dtfil,gmet,chi0,spaceComm,ec_rpa,ec_gm) 2601 2602 use m_hide_lapack, only : xginv, xheev 2603 2604 !Arguments ------------------------------------ 2605 !scalars 2606 integer,intent(in) :: iqcalc,iq,gwrpacorr,gwgmcorr,spaceComm 2607 real(dp),intent(inout) :: ec_gm 2608 type(kmesh_t),intent(in) :: Qmesh 2609 type(vcoul_t),intent(in) :: Pvc 2610 type(Datafiles_type),intent(in) :: Dtfil 2611 type(em1params_t),intent(in) :: Ep 2612 !arrays 2613 real(dp),intent(in) :: gmet(3,3) 2614 real(dp),intent(inout) :: ec_rpa(gwrpacorr) 2615 complex(gwpc),intent(inout) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega) 2616 2617 !Local variables------------------------------- 2618 !scalars 2619 integer :: ig1,ig2,ilambda,io,master,rank,nprocs,unt,ierr 2620 real(dp) :: ecorr,ecorr_gm 2621 real(dp) :: lambda 2622 logical :: q_is_gamma 2623 character(len=500) :: msg 2624 !arrays 2625 real(dp),allocatable :: z(:),zl(:),zlw(:),zw(:) 2626 complex(gwpc),allocatable :: chi0_diag(:),chitmp(:,:),chi0_diag_gm(:),chitmp_gm(:,:) 2627 real(gwpc),allocatable :: eig(:) 2628 2629 ! ************************************************************************* 2630 2631 DBG_ENTER("COLL") 2632 2633 ! initialize MPI data 2634 master=0 2635 rank = xmpi_comm_rank(spaceComm) 2636 nprocs = xmpi_comm_size(spaceComm) 2637 2638 !if (rank==master) then ! presently only master has chi0 in screening 2639 2640 ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff 2641 q_is_gamma = normv(Qmesh%ibz(:,iq),gmet,'G')<GW_TOLQ0 2642 2643 ! Calculate Gauss-Legendre quadrature knots and weights for the omega integration 2644 ABI_MALLOC(zw, (Ep%nomegaei)) 2645 ABI_MALLOC(z, (Ep%nomegaei)) 2646 call coeffs_gausslegint(zero,one,z,zw,Ep%nomegaei) 2647 2648 ! Calculate Gauss-Legendre quadrature knots and weights for the lambda integration 2649 ABI_MALLOC(zlw,(gwrpacorr)) 2650 ABI_MALLOC(zl,(gwrpacorr)) 2651 call coeffs_gausslegint(zero,one,zl,zlw,gwrpacorr) 2652 2653 2654 ABI_MALLOC(chi0_diag,(Ep%npwe)) 2655 ABI_MALLOC_OR_DIE(chitmp,(Ep%npwe,Ep%npwe), ierr) 2656 if(gwgmcorr==1) then 2657 ABI_MALLOC(chi0_diag_gm,(Ep%npwe)) 2658 ABI_MALLOC_OR_DIE(chitmp_gm,(Ep%npwe,Ep%npwe), ierr) 2659 end if 2660 2661 do io=2,Ep%nomega 2662 !if (q_is_gamma) then 2663 ! call wrtout([std_out, ab_out], "RPA: Ignoring q==0"); cycle 2664 !end if 2665 2666 if(gwrpacorr==1) then ! exact integration over the coupling constant 2667 2668 if(modulo(io-2,nprocs)/=rank) cycle ! distributing the workload 2669 2670 do ig2=1,Ep%npwe 2671 do ig1=1,Ep%npwe 2672 chitmp(ig1,ig2) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig2,iq) * chi0(ig1,ig2,io) 2673 end do !ig1 2674 end do !ig2 2675 2676 ABI_MALLOC(eig,(Ep%npwe)) 2677 call xheev('N','U',Ep%npwe,chitmp,eig) 2678 2679 do ig1=1,Ep%npwe 2680 ec_rpa(:) = ec_rpa(:) & 2681 & - zw(io-1) / ( z(io-1) * z(io-1) ) & 2682 & * Qmesh%wt(iq) * (-log( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi ) 2683 ec_gm = ec_gm & 2684 & - zw(io-1) / ( z(io-1) * z(io-1) ) & 2685 & * Qmesh%wt(iq) * ( eig(ig1) / ( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi ) 2686 end do 2687 ABI_FREE(eig) 2688 2689 else ! numerical integration over the coupling constant 2690 2691 !if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload 2692 2693 do ilambda=1,gwrpacorr 2694 if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload 2695 lambda=zl(ilambda) 2696 do ig1=1,Ep%npwe 2697 chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq)**2 * chi0(ig1,ig1,io) 2698 end do 2699 2700 if(ilambda==1 .and. gwgmcorr==1) then ! Copy v^1/2*Chi0*v^1/2 for Galitskii-Migdal 2701 chi0_diag_gm(:) = chi0_diag(:) 2702 end if 2703 2704 do ig2=1,Ep%npwe 2705 do ig1=1,Ep%npwe 2706 chitmp(ig1,ig2) = - lambda * Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io) 2707 2708 if(ilambda==1 .and. gwgmcorr==1) then ! Use lambda=1 for Galitskii-Migdal 2709 chitmp_gm(ig1,ig2) = - Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io) 2710 end if 2711 2712 end do !ig1 2713 chitmp(ig2,ig2) = chitmp(ig2,ig2) + 1.0_dp 2714 2715 if(ilambda==1 .and. gwgmcorr==1) then ! Prepare (1-v^1/2*Chi0*v^1/2) for Galitskii-Migdal 2716 chitmp_gm(ig2,ig2) = chitmp_gm(ig2,ig2) + 1.0_dp 2717 end if 2718 2719 end do !ig2 2720 call xginv(chitmp(:,:),Ep%npwe) 2721 chitmp(:,:) = matmul( chi0(:,:,io) , chitmp(:,:) ) 2722 2723 if(ilambda==1 .and. gwgmcorr==1) then ! Prepare Chi = [(1-v^1/2*Chi0*v^1/2)]^-1 * Chi0 for Galitskii-Migdal 2724 call xginv(chitmp_gm(:,:),Ep%npwe) 2725 chitmp_gm(:,:) = matmul( chi0(:,:,io) , chitmp_gm(:,:) ) 2726 end if 2727 2728 do ig1=1,Ep%npwe 2729 chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp(ig1,ig1) - chi0_diag(ig1) 2730 2731 if(ilambda==1 .and. gwgmcorr==1) then ! Prepare v^1/2*Chi*v^1/2 - v^1/2*Chi0*v^1/2 for Galitskii-Migdal 2732 chi0_diag_gm(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp_gm(ig1,ig1) - chi0_diag_gm(ig1) 2733 end if 2734 2735 end do 2736 2737 do ig1=1,Ep%npwe 2738 ec_rpa(ilambda) = ec_rpa(ilambda) & 2739 & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag(ig1) ) / (2.0_dp * pi ) 2740 2741 if(ilambda==1 .and. gwgmcorr==1) then ! Integrate [v*Chi-v*Chi0](iw) dw for Galitskii-Migdal 2742 ec_gm = ec_gm & 2743 & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag_gm(ig1) ) / (2.0_dp * pi ) 2744 end if 2745 2746 end do 2747 2748 end do ! ilambda 2749 2750 end if ! exact or numerical integration over the coupling constant 2751 2752 end do ! io 2753 2754 2755 ! Output the correlation energy when the last q-point to be calculated is reached 2756 ! This would allow for a manual parallelization over q-points 2757 if(iqcalc==Ep%nqcalc) then 2758 2759 call xmpi_sum_master(ec_rpa,master,spaceComm,ierr) 2760 call xmpi_sum_master(ec_gm,master,spaceComm,ierr) 2761 2762 if(rank==master) then 2763 ecorr = sum( zlw(:)*ec_rpa(:) ) 2764 ecorr_gm = ec_gm 2765 if (open_file(dtfil%fnameabo_rpa, msg, newunit=unt) /=0) then 2766 ABI_ERROR(msg) 2767 end if 2768 write(unt,'(a,(2x,f14.8))') '#RPA',ecorr 2769 write(msg,'(2a,(2x,f14.8))') ch10,' RPA energy [Ha] :',ecorr 2770 call wrtout([ab_out, std_out], msg) 2771 if(gwrpacorr>1) then 2772 do ilambda=1,gwrpacorr 2773 write(unt,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda) 2774 write(msg,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda) 2775 call wrtout([ab_out, std_out], msg) 2776 end do 2777 end if 2778 if(gwgmcorr==1) then ! Only exact integration over the coupling constant 2779 write(unt,'(a,(2x,f14.8))') '#GM',ecorr_gm 2780 write(msg,'(2a,(2x,f14.8))') ch10,' Galitskii-Migdal energy [Ha] :',ecorr_gm 2781 call wrtout([ab_out, std_out], msg) 2782 write(unt,'(a1)') ' ' 2783 write(msg,'(a1)') ' ' 2784 call wrtout([ab_out, std_out], msg) 2785 end if 2786 close(unt) 2787 end if 2788 2789 end if 2790 2791 if(gwgmcorr==1) then 2792 ABI_FREE(chitmp_gm) 2793 ABI_FREE(chi0_diag_gm) 2794 end if 2795 ABI_FREE(chi0_diag) 2796 ABI_FREE(chitmp) 2797 ABI_FREE(zl) 2798 ABI_FREE(zlw) 2799 ABI_FREE(z) 2800 ABI_FREE(zw) 2801 2802 DBG_EXIT("COLL") 2803 2804 end subroutine calc_rpa_functional
m_screening_driver/chi0_bksmask [ Functions ]
[ Top ] [ m_screening_driver ] [ Functions ]
NAME
chi0_bksmask
FUNCTION
Compute tables for the distribution and the storage of the wavefunctions in the SCREENING code.
INPUTS
Dtset<type(dataset_type)>=all input variables for this dataset Ep<em1params_t>=Parameters for the screening calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. nbvw = Max. number of fully/partially occupied states over spin nbcw = Max. number of unoccupied states considering the spin nprocs=Total number of MPI processors my_rank=Rank of this this processor.
OUTPUT
bks_mask(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state. keep_ur(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space. ierr=Exit status.
SOURCE
2262 subroutine chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr) 2263 2264 !Arguments ------------------------------------ 2265 !scalars 2266 integer,intent(in) :: my_rank,nprocs,nbvw,nbcw 2267 integer,intent(out) :: ierr 2268 type(Dataset_type),intent(in) :: Dtset 2269 type(em1params_t),intent(in) :: Ep 2270 type(kmesh_t),intent(in) :: Kmesh 2271 !arrays 2272 logical,intent(out) :: bks_mask(Ep%nbnds,Kmesh%nibz,Dtset%nsppol) 2273 logical,intent(out) :: keep_ur(Ep%nbnds,Kmesh%nibz,Dtset%nsppol) 2274 2275 !Local variables------------------------------- 2276 !scalars 2277 integer :: my_nspins,my_maxb,my_minb,isp,spin,nsppol,band,rank_spin,ib 2278 character(len=500) :: msg 2279 logical :: store_ur 2280 !arrays 2281 integer :: my_spins(Dtset%nsppol),nprocs_spin(Dtset%nsppol) 2282 integer,allocatable :: istart(:),istop(:) 2283 2284 ! ************************************************************************* 2285 2286 ierr=0; nsppol=Dtset%nsppol 2287 2288 my_nspins=Dtset%nsppol; my_spins= [(isp,isp=1,nsppol)] 2289 2290 ! List of spins for each node, number of processors per each spin 2291 ! and the MPI rank in the "spin" communicator. 2292 nprocs_spin = nprocs; rank_spin = my_rank 2293 2294 if (nsppol==2.and.nprocs>1) then 2295 ! Distribute spins (optimal distribution if nprocs is even) 2296 nprocs_spin(1) = nprocs/2 2297 nprocs_spin(2) = nprocs - nprocs/2 2298 my_nspins=1; my_spins(1)=1 2299 if (my_rank+1>nprocs/2) then 2300 my_spins(1)=2 2301 rank_spin = my_rank - nprocs/2 2302 end if 2303 end if 2304 2305 store_ur = (MODULO(Dtset%gwmem,10)==1) 2306 bks_mask=.FALSE.; keep_ur=.FALSE. 2307 2308 select case (Dtset%gwpara) 2309 case (1) 2310 ! Parallelization over transitions **without** memory distributions (Except for the spin). 2311 my_minb=1; my_maxb=Ep%nbnds 2312 do isp=1,my_nspins 2313 spin = my_spins(isp) 2314 bks_mask(my_minb:my_maxb,:,spin) = .TRUE. 2315 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 2316 end do 2317 2318 case (2) 2319 ! Distribute bands and spin. 2320 do isp=1,my_nspins 2321 spin = my_spins(isp) 2322 2323 if (nprocs_spin(spin) <= nbcw) then 2324 ! Distribute nbcw empty bands among nprocs_spin (block of bands without replicas). 2325 ! Bands are distributed in contiguous blocks because 2326 ! this distribution is well suited for the Hilber transform 2327 ! since each node will allocate only a smaller frequency interval 2328 ! for the spectral function whose size scales with the number of MPI nodes. 2329 ! Note it is now meaningless to distinguish gwcomp=0 or 1 since the workload is well balanced later on 2330 ABI_MALLOC(istart,(nprocs_spin(spin))) 2331 ABI_MALLOC(istop,(nprocs_spin(spin))) 2332 2333 call xmpi_split_work2_i4b(nbcw,nprocs_spin(spin),istart,istop) 2334 2335 my_minb = nbvw + istart(rank_spin+1) 2336 my_maxb = nbvw + istop (rank_spin+1) 2337 2338 ABI_FREE(istart) 2339 ABI_FREE(istop) 2340 2341 if (my_maxb - my_minb + 1 <= 0) then 2342 write(msg,'(3a,2(i0,a),2a)')& 2343 'One or more processors has zero number of bands ',ch10,& 2344 'my_minb= ',my_minb,' my_maxb= ',my_maxb,ch10,& 2345 'This is a waste, decrease the number of processors.' 2346 ABI_ERROR(msg) 2347 end if 2348 2349 bks_mask(my_minb:my_maxb,:,spin)=.TRUE. 2350 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 2351 2352 else 2353 ! New version (alternate bands with replicas if nprocs > nbcw) 2354 ! FIXME: Fix segmentation fault with Hilbert transform. 2355 do ib=1,nbcw 2356 if (xmpi_distrib_with_replicas(ib,nbcw,rank_spin,nprocs_spin(spin))) then 2357 band = ib + nbvw 2358 bks_mask(band,:,spin)=.TRUE. 2359 if (store_ur) keep_ur(band,:,spin)=.TRUE. 2360 end if 2361 end do 2362 end if 2363 2364 ! This is needed to have all the occupied states on each node. 2365 bks_mask(1:nbvw,:,spin) = .TRUE. 2366 if (store_ur) keep_ur(1:nbvw,:,spin)=.TRUE. 2367 end do ! isp 2368 2369 case default 2370 ierr = 1 2371 ABI_WARNING("Wrong value for gwpara") 2372 end select 2373 2374 end subroutine chi0_bksmask
m_screening_driver/random_stopping_power [ Functions ]
[ Top ] [ m_screening_driver ] [ Functions ]
NAME
random_stopping_power
FUNCTION
Calculate the electronic random stopping power
INPUTS
OUTPUT
SIDE EFFECTS
SOURCE
2392 subroutine random_stopping_power(iqibz,npvel,pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower) 2393 2394 use m_splines 2395 2396 !Arguments ------------------------------------ 2397 !scalars 2398 integer,intent(in) :: iqibz,npvel 2399 real(dp),intent(in) :: pvelmax(3) 2400 2401 type(em1params_t),intent(in) :: Ep 2402 type(gsphere_t),intent(in) :: Gsph_epsG0 2403 type(kmesh_t),intent(in) :: Qmesh 2404 type(vcoul_t),intent(in) :: Vcp 2405 type(crystal_t),intent(in) :: Cryst 2406 type(Datafiles_type),intent(in) :: Dtfil 2407 2408 complex(gwpc),intent(in) :: epsm1(Ep%npwe,Ep%npwe,Ep%nomega) 2409 2410 real(dp),intent(inout) :: rspower(npvel) 2411 2412 !Local variables ------------------------------ 2413 integer :: ipvel,ig 2414 integer :: iq_bz,iq_ibz,isym_q,itim_q 2415 integer :: iomega,iomegap,nomega_re 2416 integer :: unt_rsp 2417 integer,allocatable :: iomega_re(:) 2418 2419 real(dp),parameter :: zp=1.0_dp ! Hard-coded charge of the impinging particle 2420 real(dp) :: omega_p 2421 real(dp) :: im_epsm1_int(1) 2422 real(dp) :: qbz(3),qpgcart(3),qpg_red(3) 2423 real(dp) :: pvel(3,npvel),pvel_norm(npvel) 2424 real(dp) :: ypp_i(Ep%nomega) 2425 real(dp) :: vcoul(Ep%npwe) 2426 real(dp),allocatable :: im_epsm1_diag_qbz(:,:),tmp_data(:) 2427 real(dp),allocatable :: omega_re(:) 2428 2429 character(len=500) :: msg 2430 character(len=fnlen+4) :: fname 2431 2432 !************************************************************************ 2433 2434 ! 2435 ! First set up the velocities array from the input variables npvel and pvelmax(3) 2436 ! Remember pvelmax is in Cartesian coordinates and so is pvel 2437 do ipvel=1,npvel 2438 pvel(:,ipvel) = REAL(ipvel,dp) / REAL(npvel,dp) * pvelmax(:) 2439 pvel_norm(ipvel) = SQRT( SUM( pvel(:,ipvel)**2 ) ) 2440 enddo 2441 ! 2442 ! Select the purely real frequency in Ep%omega 2443 nomega_re=0 2444 do iomega=1,Ep%nomega 2445 if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then 2446 nomega_re=nomega_re+1 2447 endif 2448 enddo 2449 ABI_MALLOC(omega_re,(nomega_re)) 2450 ABI_MALLOC(iomega_re,(nomega_re)) 2451 ABI_MALLOC(im_epsm1_diag_qbz,(Ep%npwe,Ep%nomega)) 2452 ABI_MALLOC(tmp_data,(Ep%nomega)) 2453 2454 iomegap=0 2455 do iomega=1,Ep%nomega 2456 if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then 2457 iomegap=iomegap+1 2458 iomega_re(iomegap)=iomega 2459 omega_re(iomegap)=REAL(Ep%omega(iomega),dp) 2460 endif 2461 enddo 2462 2463 ! 2464 ! Loop over all the q-points in the full Brillouin zone and select only the 2465 ! ones that corresponds to the correct q-point in the irreducible wedge we are 2466 ! currently treating (index iqibz) 2467 do iq_bz=1,Qmesh%nbz 2468 2469 ! Perform the check and obtain the symmetry information 2470 call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q) 2471 if( iqibz /= iq_ibz ) cycle 2472 2473 ! Apply the symmetry operation to the diagonal of epsm1 2474 do iomega=1,nomega_re 2475 do ig=1,Ep%npwe 2476 im_epsm1_diag_qbz(Gsph_epsG0%rottb(ig,itim_q,isym_q),iomega)= AIMAG( epsm1(ig,ig,iomega_re(iomega)) ) 2477 enddo 2478 enddo 2479 ! Apply the symmetry operation to the Coulomb interaction 2480 do ig=1,Ep%npwe 2481 vcoul(Gsph_epsG0%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iqibz)**2 2482 enddo 2483 2484 ! 2485 ! Sum over G vectors 2486 do ig=1,Ep%npwe 2487 ! 2488 ! Loop over velocities 2489 do ipvel=1,npvel 2490 2491 qpg_red(:) = qbz(:) + Gsph_epsG0%gvec(:,ig) 2492 ! Transform q + G from reduced to cartesian with the symmetry operation 2493 qpgcart(:) = two_pi * Cryst%gprimd(:,1) * qpg_red(1) & 2494 & + two_pi * Cryst%gprimd(:,2) * qpg_red(2) & 2495 & + two_pi * Cryst%gprimd(:,3) * qpg_red(3) 2496 2497 ! omega_p = ( q + G ) . v 2498 omega_p = DOT_PRODUCT( qpgcart(:) , pvel(:,ipvel) ) 2499 2500 ! Check that the calculated frequency omega_p is within the omega 2501 ! range of epsm1 and thus that the interpolation will go fine 2502 if ( ABS(omega_p) > omega_re(nomega_re) ) then 2503 write(msg,'(a,e16.4,2a,e16.4)') ' freqremax is currently ',omega_re(nomega_re),ch10,& 2504 & ' increase it to at least ',omega_p 2505 ABI_WARNING(msg) 2506 endif 2507 2508 ! Perform the spline interpolation to obtain epsm1 at the desired omega = omega_p 2509 tmp_data = im_epsm1_diag_qbz(ig,:) 2510 2511 call spline( omega_re, tmp_data, nomega_re, 1.0e+32_dp, 1.0e+32_dp, ypp_i) 2512 call splint( nomega_re, omega_re, tmp_data, ypp_i, 1, (/ ABS(omega_p) /), im_epsm1_int ) 2513 2514 ! 2515 ! Apply the odd parity of Im epsm1 in omega to recover the causal response function 2516 if (omega_p<zero) im_epsm1_int(1)=-im_epsm1_int(1) 2517 2518 ! 2519 ! Calculate 4 * pi / |q+G|**2 * omega_p * Im{ epsm1_GG(q,omega_p) } 2520 ! 2521 im_epsm1_int(1) = omega_p * vcoul(ig) * im_epsm1_int(1) 2522 2523 ! Accumulate the final result without the prefactor 2524 ! (It will be included at the very end) 2525 rspower(ipvel) = rspower(ipvel) + im_epsm1_int(1) 2526 2527 end do ! end of velocity loop 2528 end do ! end G-loop 2529 2530 enddo ! end of q loop in the full BZ 2531 2532 ! If it is the last q, write down the result in the main output file and in a 2533 ! separate file _RSP (for Random Stopping Power) 2534 if (iqibz == Qmesh%nibz ) then 2535 2536 ! Multiply by the prefactors 2537 ! Note that this expression differs from Eq. (3.11) in Campillo PRB 58, 10307 (1998) [[cite:Campillo1998]]. 2538 ! A factor one half is missing in the paper. 2539 rspower(:) = - zp**2 / ( Cryst%ucvol * Qmesh%nbz * pvel_norm(:) ) * rspower(:) 2540 2541 write(msg,'(2a)') ch10,' ==== Random stopping power along Cartesian direction === ' 2542 call wrtout([ab_out, std_out], msg) 2543 write(msg,'(a,3(f12.4,2x),a)') ' ==== ',pvelmax(:),'====' 2544 call wrtout([ab_out, std_out], msg) 2545 write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) ' 2546 call wrtout([ab_out, std_out], msg) 2547 do ipvel=1,npvel 2548 write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel) 2549 call wrtout([ab_out, std_out], msg) 2550 enddo 2551 write(msg,'(2a)') ' ========================================================= ',ch10 2552 call wrtout([ab_out, std_out], msg) 2553 2554 fname=TRIM(Dtfil%filnam_ds(4))//'_RSP' 2555 if (open_file(fname,msg,newunit=unt_rsp,status='unknown',form='formatted') /= 0) then 2556 ABI_ERROR(msg) 2557 end if 2558 2559 write(msg,'(a)') '# ==== Random stopping power along Cartesian direction === ' 2560 call wrtout(unt_rsp, msg) 2561 write(msg,'(a,3(f12.4,2x))') '# ==== ',pvelmax(:) 2562 call wrtout(unt_rsp, msg) 2563 write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) ' 2564 call wrtout(unt_rsp, msg) 2565 do ipvel=1,npvel 2566 write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel) 2567 call wrtout(unt_rsp,msg) 2568 enddo 2569 close(unt_rsp) 2570 end if 2571 2572 ABI_FREE(omega_re) 2573 ABI_FREE(iomega_re) 2574 ABI_FREE(im_epsm1_diag_qbz) 2575 ABI_FREE(tmp_data) 2576 2577 end subroutine random_stopping_power
m_screening_driver/screening [ Functions ]
[ Top ] [ m_screening_driver ] [ Functions ]
NAME
screening
FUNCTION
Calculate screening and dielectric functions
INPUTS
acell(3)=length scales of primitive translations (bohr) codvsn=code version Dtfil<datafiles_type)>=variables related to file names and unit numbers. Pawang<pawang_type)>=paw angular mesh and related data Pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data Pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data Psps <type(pseudopotential_type)>=variables related to pseudopotentials Before entering the first time in screening, a significant part of Psps has been initialized: the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid, ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in the call to pspini. The next time the code enters screening, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=dimensionless real space primitive translations
OUTPUT
Output is written on the main output file. The symmetrical inverse dielectric matrix is stored in the _SCR file
SIDE EFFECTS
Dtset<type(dataset_type)>=all input variables for this dataset
NOTES
USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...Total density, potentials, ... are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
155 subroutine screening(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim) 156 157 !Arguments ------------------------------------ 158 !scalars 159 character(len=8),intent(in) :: codvsn 160 type(Datafiles_type),intent(in) :: Dtfil 161 type(Dataset_type),intent(inout) :: Dtset 162 type(Pawang_type),intent(inout) :: Pawang 163 type(Pseudopotential_type),intent(inout) :: Psps 164 !arrays 165 real(dp),intent(in) :: acell(3),rprim(3,3) 166 type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Dtset%usepaw) 167 type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Dtset%usepaw) 168 169 !Local variables ------------------------------ 170 character(len=4) :: ctype='RPA ' 171 !scalars 172 integer,parameter :: tim_fourdp4=4,NOMEGA_PRINTED=15,master=0 173 integer :: spin,ik_ibz,my_nbks 174 integer :: choice,cplex,cplex_rhoij,dim_kxcg,dim_wing,ount,omp_ncpus 175 integer :: fform_chi0,fform_em1,gnt_option,iat,ider,idir,ierr,band 176 integer :: ifft,ii,ikbz,ikxc,initialized,iomega,ios,ipert 177 integer :: iqibz,iqcalc,is_qeq0,isym,izero,ifirst,ilast 178 integer :: label,mgfftf,mgfftgw 179 integer :: nt_per_proc,work_size 180 integer :: moved_atm_inside,moved_rhor 181 integer :: nbcw,nbsc,nbvw,nkxc,nkxc1,n3xccc,optene,istep 182 integer :: nfftf,nfftf_tot,nfftgw,nfftgw_tot,ngrvdw,nhatgrdim,nprocs,nspden_rhoij 183 integer :: nscf,nzlmopt,mband 184 integer :: optcut,optgr0,optgr1,optgr2,option,approx_type,option_test,optgrad 185 integer :: optrad,optrhoij,psp_gencond,my_rank, ig 186 integer :: rhoxsp_method,comm,test_type,tordering,unt_em1,unt_susc,usexcnhat, ncerr 187 real(dp) :: compch_fft,compch_sph,domegareal,e0,ecore,ecut_eff,ecutdg_eff 188 real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp,omegaplasma,ucvol,vxcavg,gw_gsq,r_s 189 real(dp) :: alpha,rhoav,factor,ec_gm 190 real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 191 logical :: found,iscompatibleFFT,is_dfpt=.false.,use_tr,is_first_qcalc 192 logical :: add_chi0_intraband,update_energies,call_pawinit 193 character(len=10) :: string 194 character(len=500) :: msg 195 character(len=80) :: bar 196 type(ebands_t) :: ks_ebands, qp_ebands 197 type(kmesh_t) :: Kmesh,Qmesh 198 type(vcoul_t) :: Vcp 199 type(crystal_t) :: Cryst 200 type(em1params_t) :: Ep 201 type(Energies_type) :: KS_energies 202 type(gsphere_t) :: Gsph_epsG0,Gsph_wfn 203 type(Hdr_type) :: Hdr_wfk,Hdr_local 204 type(MPI_type) :: MPI_enreg_seq 205 type(Pawfgr_type) :: Pawfgr 206 type(hscr_t) :: Hem1,Hchi0 207 type(wfdgw_t) :: Wfd,Wfdf 208 type(spectra_t) :: spectra 209 type(chi_t) :: chihw 210 type(wvl_data) :: wvl_dummy 211 character(len=nctk_slen) :: wing_shape 212 !arrays 213 integer :: ibocc(Dtset%nsppol),ngfft_gw(18),ngfftc(18),ngfftf(18) 214 integer,allocatable :: irottb(:,:),ktabr(:,:),ktabrf(:,:),l_size_atm(:) 215 integer,allocatable :: ks_vbik(:,:),ks_occ_idx(:,:),qp_vbik(:,:),nband(:,:) 216 integer,allocatable :: nq_spl(:),nlmn_atm(:),gw_gfft(:,:) 217 real(dp) :: gmet(3,3),gprimd(3,3),k0(3),qtmp(3),rmet(3,3),rprimd(3,3),tsec(2),strsxc(6) 218 real(dp),allocatable :: igwene(:,:,:),chi0_sumrule(:),ec_rpa(:),rspower(:) 219 real(dp),allocatable :: nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:) 220 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhor_p(:,:),rhor_kernel(:,:),taur(:,:) 221 real(dp),allocatable :: z(:),zw(:),grchempottn(:,:),grewtn(:,:),grvdw(:,:),kxc(:,:),qmax(:) 222 real(dp),allocatable :: ks_vhartr(:),vpsp(:),ks_vtrial(:,:),ks_vxc(:,:),xccc3d(:) 223 complex(gwpc),allocatable :: arr_99(:,:),kxcg(:,:),fxc_ADA(:,:,:) 224 complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:) 225 complex(dpc),allocatable :: chi0_head(:,:,:), chi0_lwing(:,:,:), chi0_uwing(:,:,:) 226 real(dp),allocatable :: rwork_wing(:,:,:,:) 227 complex(dpc),allocatable :: chi0intra_lwing(:,:,:),chi0intra_uwing(:,:,:),chi0intra_head(:,:,:) 228 complex(gwpc),allocatable,target :: chi0(:,:,:),chi0intra(:,:,:) 229 complex(gwpc),ABI_CONTIGUOUS pointer :: epsm1(:,:,:) 230 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 231 character(len=80) :: title(2) 232 character(len=fnlen) :: gw_fname,wfk_fname,lwl_fname 233 type(littlegroup_t),pointer :: Ltg_q(:) 234 type(Paw_an_type),allocatable :: Paw_an(:) 235 type(Paw_ij_type),allocatable :: Paw_ij(:) 236 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 237 type(Pawrhoij_type),allocatable :: Pawrhoij(:),prev_Pawrhoij(:) 238 type(pawpwff_t),allocatable :: Paw_pwff(:) 239 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 240 type(plowannier_type) :: wanbz,wanibz,wanibz_in 241 242 #ifdef __HAVE_GREENX 243 integer :: gap_err 244 real(dp) :: te_min, te_max 245 type(gaps_t) :: gaps 246 real(dp),allocatable :: tau_mesh(:), tau_wgs(:), iw_mesh(:), iw_wgs(:) 247 real(dp),allocatable :: t2w_cos_wgs(:,:), w2t_cos_wgs(:,:), t2w_sin_wgs(:,:) 248 real(dp) :: ft_max_error(3), cosft_duality_error 249 #endif 250 251 !************************************************************************ 252 253 DBG_ENTER("COLL") 254 255 call timab(301,1,tsec) ! overall time 256 call timab(302,1,tsec) ! screening(init 257 258 write(msg,'(6a)')& 259 ' SCREENING: Calculation of the susceptibility and dielectric matrices ',ch10,ch10,& 260 ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,& 261 ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.' 262 call wrtout([ab_out, std_out], msg) 263 264 if(dtset%ucrpa>0) then 265 write(msg,'(6a)')ch10,& 266 ' cRPA Calculation: The calculation of the polarisability is constrained (ucrpa/=0)',ch10 267 call wrtout([ab_out, std_out], msg) 268 end if 269 #if defined HAVE_GW_DPC 270 if (gwpc/=8) then 271 write(msg,'(6a)')ch10,& 272 ' Number of bytes for double precision complex /=8 ',ch10,& 273 ' Cannot continue due to kind mismatch in BLAS library ',ch10,& 274 ' Some BLAS interfaces are not generated by abilint ' 275 ABI_ERROR(msg) 276 end if 277 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 278 #else 279 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 280 #endif 281 call wrtout([ab_out, std_out], msg) 282 283 ! === Initialize MPI variables, and parallelization level === 284 ! gwpara: 0--> sequential run, 1--> parallelism over k-points, 2--> parallelism over bands. 285 ! gwpara==2, each node has both fully and partially occupied states while conduction bands are divided 286 comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 287 288 if (my_rank == master) then 289 wfk_fname = dtfil%fnamewffk 290 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 291 ABI_ERROR(msg) 292 end if 293 end if 294 call xmpi_bcast(wfk_fname, master, comm, ierr) 295 296 ! Some variables need to be initialized/nullify at start 297 call energies_init(KS_energies) 298 usexcnhat=0 299 300 call mkrdim(acell,rprim,rprimd) 301 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 302 303 !=== Define FFT grid(s) sizes === 304 ! Be careful! This mesh is only used for densities and potentials. It is NOT the (usually coarser) 305 ! GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 306 ! See also NOTES in the comments at the beginning of this file. 307 ! NOTE: The mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 308 309 k0(:)=zero 310 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 311 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 312 313 call print_ngfft(ngfftf,'Dense FFT mesh used for densities and potentials') 314 nfftf_tot=PRODUCT(ngfftf(1:3)) 315 316 ! We can intialize MPI_enreg and fft distrib here, now ngfft are known 317 call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for the sequential part. 318 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 319 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 320 321 !============================================= 322 !==== Open and read pseudopotential files ==== 323 !============================================= 324 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 325 326 ! === Initialize dimensions and basic objects === 327 call setup_screening(codvsn,acell,rprim,wfk_fname,Dtset,Psps,Pawtab,& 328 ngfft_gw,Hdr_wfk,Hdr_local,Cryst,Kmesh,Qmesh,ks_ebands,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm) 329 330 call timab(302,2,tsec) ! screening(init) 331 call print_ngfft(ngfft_gw,'FFT mesh used for oscillator strengths') 332 333 nfftgw_tot=PRODUCT(ngfft_gw(1:3)) 334 mgfftgw =MAXVAL (ngfft_gw(1:3)) 335 nfftgw =nfftgw_tot ! no FFT // 336 337 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 338 KS_energies%e_corepsp=ecore/Cryst%ucvol 339 340 !========================== 341 !=== PAW initialization === 342 !========================== 343 if (Dtset%usepaw==1) then 344 call timab(315,1,tsec) ! screening(pawin 345 346 call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,Cryst%xred) 347 348 ABI_MALLOC(Pawrhoij,(Cryst%natom)) 349 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 350 & nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc) 351 call pawrhoij_alloc(Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,& 352 & Cryst%typat,pawtab=Pawtab) 353 354 ! Initialize values for several basic arrays stored in Pawinit 355 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 356 357 ! Test if we have to call pawinit 358 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 359 360 if (psp_gencond==1.or.call_pawinit) then 361 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 362 call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,& 363 & Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,& 364 & Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%ixc,Dtset%usepotzero) 365 366 ! Update internal values 367 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 368 else 369 if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2 370 if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2 371 end if 372 Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore) 373 374 ! Initialize optional flags in Pawtab to zero 375 ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed) 376 Pawtab(:)%has_nabla = 0 377 Pawtab(:)%usepawu = 0 378 Pawtab(:)%useexexch = 0 379 Pawtab(:)%exchmix =zero 380 Pawtab(:)%lamb_shielding = zero 381 382 ! Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit. 383 ! TODO solve problem with memory leak and clean this part as well as the associated flag 384 call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab) 385 386 call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,rprimd,Cryst%symrec,Pawang%zarot) 387 388 ! Initialize and compute data for DFT+U. 389 ! paw_dmft%use_dmft=dtset%usedmft 390 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 391 & is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,& 392 & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa) 393 394 if (my_rank == master) call pawtab_print(Pawtab) 395 396 ! Get Pawrhoij from the header of the WFK file. 397 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 398 399 ! Re-symmetrize rhoij. 400 ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly 401 choice=1; optrhoij=1; ipert=0; idir=0 402 !call pawrhoij_symrhoij(Pawrhoij,Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,& 403 !& Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,& 404 !& Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 405 ! 406 ! Evaluate form factors for the radial part of phi.phj-tphi.tphj === 407 ! rhoxsp_method=1 ! Arnaud-Alouani 408 ! rhoxsp_method=2 ! Shiskin-Kresse 409 rhoxsp_method=2 410 411 ! At least for ucrpa, the Arnaud Alouani is always a better choice but needs a larger cutoff 412 if(dtset%ucrpa>0) rhoxsp_method=1 413 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 414 415 ABI_MALLOC(gw_gfft,(3,nfftgw_tot)) 416 call get_gftt(ngfft_gw,(/zero,zero,zero/),gmet,gw_gsq,gw_gfft) 417 ABI_FREE(gw_gfft) 418 419 ! Set up q grids, make qmax 20% larger than largest expected: 420 ABI_MALLOC(nq_spl,(Psps%ntypat)) 421 ABI_MALLOC(qmax,(Psps%ntypat)) 422 nq_spl = Psps%mqgrid_ff 423 qmax = SQRT(gw_gsq)*1.2d0 !qmax = Psps%qgrid_ff(Psps%mqgrid_ff) 424 ABI_MALLOC(Paw_pwff,(Psps%ntypat)) 425 426 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 427 ABI_FREE(nq_spl) 428 ABI_FREE(qmax) 429 430 ! Variables/arrays related to the fine FFT grid 431 ABI_MALLOC(nhat,(nfftf,Dtset%nspden)) 432 nhat=zero; cplex=1 433 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 434 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 435 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat) 436 ABI_FREE(l_size_atm) 437 compch_fft=greatest_real 438 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 439 ! * 0 --> Vloc in atomic data is Vbare (Blochl s formulation) 440 ! * 1 --> Vloc in atomic data is VH(tnzc) (Kresse s formulation) 441 write(msg,'(a,i3)')' screening : using usexcnhat = ',usexcnhat 442 call wrtout(std_out, msg) 443 444 ! Identify parts of the rectangular grid where the density has to be calculated. 445 optcut=0;optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 446 if (Dtset%pawcross==1) optrad=1 447 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 448 449 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 450 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 451 452 call timab(315,2,tsec) ! screening(pawin 453 else 454 ! allocate empty structure for the sake of -fcheck=all... 455 ABI_MALLOC(Paw_pwff,(0)) 456 ABI_MALLOC(Pawrhoij,(0)) 457 ABI_MALLOC(Pawfgrtab,(0)) 458 end if ! End of PAW initialization. 459 460 ! Consistency check and additional stuff done only for GW with PAW. 461 ABI_MALLOC(Paw_onsite,(Cryst%natom)) 462 if (Dtset%usepaw==1) then 463 if (Dtset%ecutwfn < Dtset%ecut) then 464 write(msg,"(5a)")& 465 "WARNING - ",ch10,& 466 " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 467 " an excessive truncation of the planewave basis set can lead to unphysical results." 468 call wrtout(ab_out, msg) 469 end if 470 471 ABI_CHECK(Dtset%useexexch==0,"LEXX not yet implemented in GW") 472 ABI_CHECK(Dtset%usedmft==0,"DMFT + GW not available") 473 474 if (Dtset%pawcross==1) then 475 optgrad=1 476 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,& 477 Cryst%xcart,Pawtab,Pawrad,Pawfgrtab,optgrad) 478 end if 479 end if 480 481 ! Allocate these arrays anyway, since they are passed to subroutines. 482 if (.not.allocated(nhat)) then 483 ABI_MALLOC(nhat,(nfftf,0)) 484 end if 485 486 call timab(316,1,tsec) ! screening(wfs 487 488 !===================================================== 489 !=== Prepare the distribution of the wavefunctions === 490 !===================================================== 491 ! valence and partially occupied are replicate on each node while conduction bands are MPI distributed. 492 ! This method is mandatory if gwpara==2 and/or we are using awtr==1 or the spectral method. 493 ! If awtr==1, we evaluate chi0 taking advantage of time-reversal (speed-up~2) 494 ! Useful indices: 495 ! nbvw = Max. number of fully/partially occupied states over spin 496 ! nbcw = Max. number of unoccupied states considering the spin 497 !TODO: 498 ! Here for semiconducting systems we have to be sure that each processor has all the 499 ! states considered in the SCGW, moreover nbsc<nbvw 500 ! in case of SCGW vale and conduction has to be recalculated to avoid errors 501 ! if a metal becomes semiconductor or viceversa. 502 ! Ideally nbvw should include only the states v such that the transition 503 ! c-->v is taken into account in cchi0 (see GW_TOLDOCC). In the present implementation 504 505 ABI_MALLOC(ks_occ_idx,(ks_ebands%nkpt, ks_ebands%nsppol)) 506 ABI_MALLOC(ks_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol)) 507 ABI_MALLOC(qp_vbik ,(ks_ebands%nkpt, ks_ebands%nsppol)) 508 509 call ebands_update_occ(ks_ebands, Dtset%spinmagntarget, prtvol=0) 510 ks_occ_idx = ebands_get_occupied(ks_ebands,tol8) ! tol8 to be consistent when the density 511 ks_vbik = ebands_get_valence_idx(ks_ebands) 512 513 ibocc(:)=MAXVAL(ks_occ_idx(:,:),DIM=1) ! Max occupied band index for each spin. 514 ABI_FREE(ks_occ_idx) 515 516 use_tr=.FALSE.; nbvw=0 517 if (Dtset%gwpara==2.or.Ep%awtr==1.or.Dtset%spmeth>0) then 518 use_tr=.TRUE. 519 nbvw=MAXVAL(ibocc) 520 nbcw=Ep%nbnds-nbvw 521 write(msg,'(4a,i0,2a,i0,2a,i0,a)')ch10,& 522 '- screening: taking advantage of time-reversal symmetry ',ch10,& 523 '- Maximum band index for partially occupied states nbvw = ',nbvw,ch10,& 524 '- Remaining bands to be divided among processors nbcw = ',nbcw,ch10,& 525 '- Number of bands treated by each node ~',nbcw/nprocs,ch10 526 call wrtout(ab_out, msg) 527 if (Cryst%timrev/=2) then 528 ABI_ERROR('Time-reversal cannot be used since cryst%timrev/=2') 529 end if 530 end if 531 532 mband=Ep%nbnds 533 ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol)) 534 nband=mband 535 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol)) 536 ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol)) 537 bks_mask=.FALSE.; keep_ur=.FALSE. 538 539 ! autoparal section 540 if (dtset%max_ncpus /=0) then 541 ount = ab_out 542 ! Temporary table needed to estimate memory 543 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 544 if (Dtset%usepaw==1) then 545 do iat=1,Cryst%natom 546 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 547 end do 548 end if 549 550 write(ount,'(a)')"--- !Autoparal" 551 write(ount,"(a)")"# Autoparal section for Screening runs" 552 write(ount,"(a)") "info:" 553 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 554 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 555 write(ount,"(a,i0)")" gwpara: ",dtset%gwpara 556 write(ount,"(a,i0)")" nkpt: ",dtset%nkpt 557 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 558 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 559 write(ount,"(a,i0)")" nbnds: ",Ep%nbnds 560 561 work_size = nbvw * nbcw * Kmesh%nibz**2 * Dtset%nsppol 562 563 ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI. 564 nonscal_mem = (two*gwpc*Ep%npwe**2*(Ep%nomega*b2Mb)) * 1.1_dp 565 566 ! List of configurations. 567 ! Assuming an OpenMP implementation with perfect speedup! 568 write(ount,"(a)")"configurations:" 569 570 do ii=1,dtset%max_ncpus 571 nt_per_proc = 0 572 eff = HUGE(one) 573 max_wfsmem_mb = zero 574 575 do my_rank=0,ii-1 576 call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,ii,bks_mask,keep_ur,ierr) 577 if (ierr /= 0) exit 578 nt_per_proc = MAX(nt_per_proc, COUNT(bks_mask(1:nbvw,:,:)) * COUNT(bks_mask(nbvw+1:,:,:))) 579 eff = MIN(eff, (one * work_size) / (ii * nt_per_proc)) 580 581 ! Memory needed for Fourier components ug. 582 my_nbks = COUNT(bks_mask) 583 ug_mem = two*gwpc*Dtset%nspinor*Ep%npwwfn*my_nbks*b2Mb 584 585 ! Memory needed for real space ur. 586 ur_mem = two*gwpc*Dtset%nspinor*nfftgw*COUNT(keep_ur)*b2Mb 587 588 ! Memory needed for PAW projections Cprj 589 cprj_mem = zero 590 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 591 592 max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem) 593 end do 594 if (ierr /= 0) cycle 595 596 ! Add the non-scalable part and increase by 10% to account for other datastructures. 597 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 598 599 do omp_ncpus=1,xomp_get_max_threads() 600 write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 601 write(ount,"(a,i0)")" mpi_ncpus: ",ii 602 write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 603 write(ount,"(a,f12.9)")" efficiency: ",eff 604 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 605 end do 606 end do 607 write(ount,'(a)')"..." 608 609 ABI_FREE(nlmn_atm) 610 ABI_ERROR_NODUMP("aborting now") 611 else 612 call chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr) 613 end if 614 615 ! Initialize the wf descriptor (allocate %ug and %ur if required). 616 617 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,& 618 Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,& 619 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 620 621 if (Dtset%pawcross==1) then 622 call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,mband,nband,Ep%nkibz,Dtset%nsppol,bks_mask,& 623 Dtset%nspden,Dtset%nspinor,dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_gw,& 624 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 625 end if 626 627 ABI_FREE(bks_mask) 628 ABI_FREE(nband) 629 ABI_FREE(keep_ur) 630 631 call wfd%print(mode_paral='PERS') 632 !FIXME: Rewrite the treatment of use_tr branches in cchi0 ... 633 !Use a different nbvw for each spin. 634 !Now use_tr means that one can use time-reversal symmetry. 635 636 !================================================== 637 !==== Read KS band structure from the KSS file ==== 638 !================================================== 639 call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname)) 640 641 if (Dtset%pawcross==1) then 642 call wfdgw_copy(Wfd, Wfdf) 643 call wfdf%change_ngfft(Cryst,Psps,ngfftf) 644 end if 645 646 ! This test has been disabled (too expensive!) 647 if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 648 649 call timab(316,2,tsec) ! screening(wfs 650 call timab(319,1,tsec) ! screening(1) 651 652 if (Cryst%nsym/=Dtset%nsym .and. Dtset%usepaw==1) then 653 ABI_ERROR('Cryst%nsym/=Dtset%nsym, check pawinit and pawrhoij_symrhoij') 654 end if 655 656 ! Get the FFT index of $ (R^{-1}(r-\tau)) $ 657 ! S= $\transpose R^{-1}$ and k_BZ = S k_IBZ 658 ! irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk. 659 ABI_MALLOC(irottb, (nfftgw, Cryst%nsym)) 660 call rotate_FFT_mesh(Cryst%nsym, Cryst%symrel, Cryst%tnons, ngfft_gw, irottb, iscompatibleFFT) 661 662 ABI_MALLOC(ktabr,(nfftgw,Kmesh%nbz)) 663 do ikbz=1,Kmesh%nbz 664 isym=Kmesh%tabo(ikbz) 665 do ifft=1,nfftgw 666 ktabr(ifft,ikbz)=irottb(ifft,isym) 667 end do 668 end do 669 ABI_FREE(irottb) 670 671 if (Dtset%usepaw==1 .and. Dtset%pawcross==1) then 672 ABI_MALLOC(irottb,(nfftf,Cryst%nsym)) 673 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT) 674 675 ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz)) 676 do ikbz=1,Kmesh%nbz 677 isym=Kmesh%tabo(ikbz) 678 do ifft=1,nfftf 679 ktabrf(ifft,ikbz)=irottb(ifft,isym) 680 end do 681 end do 682 ABI_FREE(irottb) 683 else 684 ABI_MALLOC(ktabrf,(0,0)) 685 end if 686 687 !=== Compute structure factor phases and large sphere cut-off === 688 !WARNING cannot use Dtset%mgfft, this has to be checked better 689 !mgfft=MAXVAL(ngfftc(:)) 690 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom)) 691 !write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf 692 !if (Dtset%mgfftdg/=mgfftf) write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf" 693 ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom)) 694 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom)) 695 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 696 697 if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then 698 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 699 else 700 ph1df(:,:)=ph1d(:,:) 701 end if 702 703 ! Initialize qp_ebands using KS bands 704 ! In case of SCGW, update qp_ebands using the QPS file. 705 call ebands_copy(ks_ebands, qp_ebands) 706 707 call timab(319,2,tsec) ! screening(1) 708 709 !============================ 710 !==== Self-consistent GW ==== 711 !============================ 712 if (Ep%gwcalctyp>=10) then 713 call timab(304,1,tsec) ! KS => QP; [wfrg] 714 715 ! Initialize with KS eigenvalues and eigenfunctions. 716 ABI_MALLOC(m_ks_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 717 m_ks_to_qp = czero 718 do spin=1,Wfd%nsppol 719 do ik_ibz=1,Wfd%nkibz 720 do band=1,Wfd%nband(ik_ibz,spin) 721 m_ks_to_qp(band,band,:,:) = cone 722 end do 723 end do 724 end do 725 726 ! Read unitary transformation and QP energies. 727 ! TODO switch on the renormalization of n in screening, QPS should report bdgw 728 ABI_MALLOC(rhor_p,(nfftf,Dtset%nspden)) 729 ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Psps%usepaw)) 730 731 call rdqps(qp_ebands,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,& 732 nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,rhor_p,prev_Pawrhoij) 733 734 ABI_FREE(rhor_p) 735 ABI_FREE(prev_Pawrhoij) 736 737 ! FIXME this is to preserve the old implementation for the head and the wings in ccchi0q0 738 ! But has to be rationalized 739 if (dtset%use_oldchi == 1) then 740 ks_ebands%eig = qp_ebands%eig 741 end if 742 743 ! Calculate new occ. factors and fermi level. 744 call ebands_update_occ(qp_ebands, Dtset%spinmagntarget) 745 qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands) 746 747 ! === Update only the wfg treated with GW === 748 ! For PAW update and re-symmetrize cprj in the full BZ, TODO add rotation in spinor space 749 if (nscf/=0) call wfd%rotate(Cryst,m_ks_to_qp) 750 751 ABI_FREE(m_ks_to_qp) 752 call timab(304,2,tsec) 753 end if ! gwcalctyp>=10 754 755 call timab(305,1,tsec) ! screening(densit 756 ! 757 !=== In case update the eigenvalues === 758 !* Either use a scissor operator or an external GW file. 759 gw_fname = "__in.gw__" 760 update_energies = file_exists(gw_fname) 761 762 if (ABS(Ep%mbpt_sciss)>tol6) then 763 write(msg,'(5a,f7.3,a)')& 764 ' screening : performing a first self-consistency',ch10,& 765 ' update of the energies in W by a scissor operator',ch10,& 766 ' applying a scissor operator of [eV] : ',Ep%mbpt_sciss*Ha_eV,ch10 767 call wrtout([ab_out, std_out], msg) 768 call ebands_apply_scissors(qp_ebands,Ep%mbpt_sciss) 769 else if (update_energies) then 770 write(msg,'(4a)')& 771 ' screening : performing a first self-consistency',ch10,& 772 ' update of the energies in W by a previous GW calculation via GW file: ',TRIM(gw_fname) 773 call wrtout([ab_out, std_out], msg) 774 ABI_MALLOC(igwene,(qp_ebands%mband, qp_ebands%nkpt, qp_ebands%nsppol)) 775 call rdgw(qp_ebands,gw_fname,igwene,extrapolate=.FALSE.) 776 !call rdgw(qp_ebands,gw_fname,igwene,extrapolate=.TRUE.) 777 ABI_FREE(igwene) 778 call ebands_update_occ(qp_ebands, Dtset%spinmagntarget) 779 end if 780 781 !======================== 782 !=== COMPUTE DENSITY ==== 783 !======================== 784 !* Evaluate PW part (complete charge in case of NC pseudos) 785 !TODO this part has to be rewritten. If I decrease the tol on the occupations 786 !I have to code some MPI stuff also if use_tr==.TRUE. 787 788 ABI_MALLOC(rhor,(nfftf,Dtset%nspden)) 789 ABI_MALLOC(taur,(nfftf,Dtset%nspden*Dtset%usekden)) 790 791 call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, rhor) 792 if (Dtset%usekden==1) call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, taur, optcalc=1) 793 794 call timab(305,2,tsec) ! screening(densit 795 796 nhatgrdim = 0 797 if (Dtset%usepaw==1) then ! Additional computation for PAW. 798 call timab(320,1,tsec) ! screening(paw 799 800 ! Add the compensation charge to the PW density. 801 nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 802 cplex=1; ider=2*nhatgrdim; izero=0 803 if (nhatgrdim>0) then 804 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,3)) 805 else 806 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0)) 807 end if 808 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,Cryst%gprimd,& 809 & Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 810 & Pawfgrtab,nhatgr,nhat,Pawrhoij,Pawrhoij,Pawtab,k0,Cryst%rprimd,Cryst%ucvol,dtset%usewvl,Cryst%xred) 811 812 ! === Evaluate onsite energies, potentials, densities === 813 ! * Initialize variables/arrays related to the PAW spheres. 814 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 815 cplex=1 816 ABI_MALLOC(Paw_ij,(Cryst%natom)) 817 call paw_ij_nullify(Paw_ij) 818 call paw_ij_init(Paw_ij,cplex,Dtset%nspinor,Wfd%nsppol,& 819 & Wfd%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 820 & has_dij=1,has_dijhartree=1,has_exexch_pot=1,has_pawu_occ=1) 821 822 nkxc1=0 823 ABI_MALLOC(Paw_an,(Cryst%natom)) 824 call paw_an_nullify(Paw_an) 825 call paw_an_init(Paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 826 & cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0) 827 828 nzlmopt=-1; option=0; compch_sph=greatest_real 829 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert,Dtset%ixc,& 830 & Cryst%natom,Cryst%natom,Dtset%nspden,Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,Paw_an,Paw_an,& 831 & Paw_ij,Pawang,Dtset%pawprtvol,Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,& 832 & Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,Cryst%ucvol,Psps%znuclpsp) 833 call timab(320,2,tsec) ! screening(paw 834 else 835 ABI_MALLOC(Paw_ij,(0)) 836 ABI_MALLOC(Paw_an,(0)) 837 end if ! usepaw 838 839 call timab(321,1,tsec) ! screening(2) 840 841 !JB : Should be remove : cf. l 839 842 if (.not.allocated(nhatgr)) then 843 ABI_MALLOC(nhatgr,(nfftf,Dtset%nspden,0)) 844 end if 845 846 call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,rhor,ucvol,& 847 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,omegaplasma) 848 849 !For PAW, add the compensation charge the FFT mesh, then get rho(G). 850 if (Dtset%usepaw==1) rhor(:,:)=rhor(:,:)+nhat(:,:) 851 852 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,rhor,ucvol=ucvol) 853 if(Dtset%usekden==1)then 854 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,taur,ucvol=ucvol,optrhor=1) 855 end if 856 857 if (dtset%gwgamma>0 .or. dtset%gwgamma==-11) then 858 ABI_MALLOC(rhor_kernel,(nfftf,Dtset%nspden)) 859 end if 860 861 ABI_MALLOC(rhog,(2,nfftf)) 862 call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp4) 863 864 !The following steps have been gathered in the setvtr routine: 865 !- get Ewald energy and Ewald forces 866 !- compute local ionic pseudopotential vpsp 867 !- eventually compute 3D core electron density xccc3d 868 !- eventually compute vxc and vhartr 869 !- set up ks_vtrial 870 !************************************************************** 871 !**** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 872 !************************************************************** 873 874 ngrvdw=0 875 ABI_MALLOC(grvdw,(3,ngrvdw)) 876 ABI_MALLOC(grchempottn,(3,Cryst%natom)) 877 ABI_MALLOC(grewtn,(3,Cryst%natom)) 878 nkxc=0 879 if (Dtset%nspden==1) nkxc=2 880 if (Dtset%nspden>=2) nkxc=3 ! check GGA and spinor that is messy !!! 881 ! If MGGA, fxc and kxc are not available and we dont need them for the screening part (for now ...) 882 if (Dtset%ixc<0 .and. libxc_functionals_ismgga()) nkxc=0 883 if (nkxc/=0) then 884 ABI_MALLOC(kxc,(nfftf,nkxc)) 885 end if 886 887 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 888 ABI_MALLOC(xccc3d,(n3xccc)) 889 ABI_MALLOC(ks_vhartr,(nfftf)) 890 ABI_MALLOC(ks_vtrial,(nfftf,Dtset%nspden)) 891 ABI_MALLOC(vpsp,(nfftf)) 892 ABI_MALLOC(ks_vxc,(nfftf,Dtset%nspden)) 893 894 optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1 895 call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn, & 896 & grewtn,grvdw,gsqcutf_eff,istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq, & 897 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,Cryst%ntypat,& 898 & Psps%n1xccc,n3xccc,optene,Pawang,Pawrad,Pawrhoij,Pawtab,ph1df,Psps,rhog,rhor, & 899 & Cryst%rmet,Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc, & 900 & vxcavg,wvl_dummy,xccc3d,Cryst%xred,taur=taur) 901 902 if (nkxc/=0) then 903 ABI_FREE(kxc) 904 end if 905 ABI_FREE(grchempottn) 906 ABI_FREE(grewtn) 907 ABI_FREE(grvdw) 908 ABI_FREE(xccc3d) 909 910 !============================ 911 !==== Compute KS PAW Dij ==== 912 !============================ 913 if (Dtset%usepaw==1) then 914 call timab(561,1,tsec) 915 916 ! Calculate unsymmetrized Dij. 917 cplex=1; ipert=0; idir=0 918 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert,& 919 & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 920 & Dtset%nspden,Cryst%ntypat,Paw_an,Paw_ij,Pawang,Pawfgrtab,Dtset%pawprtvol,& 921 & Pawrad,Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,k0,Dtset%spnorbscl,& 922 & Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,& 923 & nucdipmom=Dtset%nucdipmom) 924 925 ! Symmetrize KS Dij 926 #if 0 927 call symdij(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,& 928 & Cryst%nsym,Cryst%ntypat,0,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,& 929 & Cryst%symrec) 930 #else 931 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert,Cryst%natom,Cryst%natom,& 932 & Cryst%nsym,Cryst%ntypat,Paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,& 933 & Cryst%symrec) 934 #endif 935 ! 936 ! Output of the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 937 call pawprt(Dtset,Cryst%natom,Paw_ij,Pawrhoij,Pawtab) 938 call timab(561,2,tsec) 939 end if 940 ! 941 ! Calculate frequency mesh. 942 ! First omega is always zero without broadening. 943 ! FIXME what about metals? I think we should add eta, 944 ! this means we need to know if the system is metallic, for example using occopt 945 ! MS Modified to account for non-zero starting frequency (19-11-2010) 946 ! MS Modified for tangent grid (07-01-2011) 947 ABI_MALLOC(Ep%omega, (Ep%nomega)) 948 Ep%omega(1) = CMPLX(Ep%omegaermin, zero, kind=dpc) 949 !Ep%omega(1) = j_dpc * 4.372035E-01 * eV_Ha 950 !Ep%omega(1) = j_dpc * 4.372035E-01 * eV_Ha * 2 951 952 if (Ep%nomegaer > 1) then 953 ! Avoid division by zero. 954 if (Dtset%gw_frqre_tangrid == 0 .and. Dtset%gw_frqre_inzgrid == 0) then 955 domegareal = (Ep%omegaermax -Ep%omegaermin) / (Ep%nomegaer -1) 956 do iomega=2,Ep%nomegaer 957 Ep%omega(iomega)=CMPLX(Ep%omegaermin+(iomega-1)*domegareal,zero,kind=dpc) 958 end do 959 else if (Dtset%gw_frqre_tangrid == 1.and. Dtset%gw_frqre_inzgrid == 0) then 960 ! We have tangent transformed grid 961 ABI_WARNING('EXPERIMENTAL - Using tangent transform grid for contour deformation.') 962 Ep%omegaermax = Dtset%cd_max_freq 963 Ep%omegaermin = zero 964 ifirst=1; ilast=Ep%nomegaer 965 if (Dtset%cd_subset_freq(1)/=0) then ! Only a subset of frequencies is being calculated 966 ifirst=Dtset%cd_subset_freq(1); ilast=Dtset%cd_subset_freq(2) 967 end if 968 factor = Dtset%cd_halfway_freq/TAN(pi*quarter) 969 ! Important: here nfreqre is used because the step is set by the original grid 970 domegareal=(ATAN(Ep%omegaermax/factor)*two*piinv)/(Dtset%nfreqre-1) ! Stepsize in transformed variable 971 do iomega=1,Ep%nomegaer 972 Ep%omega(iomega)=CMPLX(factor*TAN((iomega+ifirst-2)*domegareal*pi*half),zero,kind=dpc) 973 end do 974 Ep%omegaermin = REAL(Ep%omega(1)) 975 Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer)) 976 else if (Dtset%gw_frqre_tangrid==0.and.Dtset%gw_frqre_inzgrid==1) then 977 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 978 domegareal=one/(Ep%nomegaer) 979 do iomega=1,Ep%nomegaer 980 factor = (iomega-1)*domegareal 981 Ep%omega(iomega)=CMPLX(e0*factor/(one-factor),zero,kind=dpc) 982 end do 983 Ep%omegaermin = REAL(Ep%omega(1)) 984 Ep%omegaermax = REAL(Ep%omega(Ep%nomegaer)) 985 else 986 ABI_ERROR('Error in specification of real frequency grid') 987 end if 988 end if 989 990 if (Ep%plasmon_pole_model .and. Ep%nomega == 2) then 991 e0= Dtset%ppmfrq; if (e0 < 0.1d-4) e0 = omegaplasma 992 Ep%omega(2)=CMPLX(zero,e0,kind=dpc) 993 end if 994 995 ! For AC, use Gauss-Legendre quadrature method. 996 ! - Replace $ \int_0^\infty dx f(x) $ with $ \int_0^1 dz f(1/z - 1)/z^2 $. 997 ! - Note that the grid is not log as required by CD thus we cannot use the same SCR file. 998 if (Ep%analytic_continuation) then 999 ABI_MALLOC(z, (Ep%nomegaei)) 1000 ABI_MALLOC(zw, (Ep%nomegaei)) 1001 call coeffs_gausslegint(zero, one, z, zw, Ep%nomegaei) 1002 do iomega=1,Ep%nomegaei 1003 Ep%omega(Ep%nomegaer + iomega) = CMPLX(zero, one/z(iomega) - one, kind=dpc) 1004 end do 1005 ABI_FREE(z) 1006 ABI_FREE(zw) 1007 1008 if (dtset%userie == 4242) then 1009 call wrtout(std_out, sjoin("userie == 4242 --> Using minimax mesh with ntau:", itoa(dtset%nfreqim))) 1010 gaps = ebands_get_gaps(ks_ebands, gap_err) 1011 ABI_CHECK(gap_err == 0, "gap_err") 1012 ! ================================ 1013 ! Setup tau/omega mesh and weights 1014 ! ================================ 1015 ! Compute min/max transition energy taking into account nsppol if any. 1016 te_min = minval(gaps%cb_min - gaps%vb_max) 1017 te_max = maxval(ks_ebands%eig(mband,:,:) - ks_ebands%eig(1,:,:)) 1018 if (te_min <= tol6) then 1019 te_min = tol6 1020 ABI_ERROR("System is metallic or with a very small fundamental gap!") 1021 end if 1022 1023 call gx_minimax_grid(dtset%nfreqim, te_min, te_max, & ! in 1024 tau_mesh, tau_wgs, & ! all these args are out and allocated by the routine. 1025 iw_mesh, iw_wgs, & 1026 t2w_cos_wgs, w2t_cos_wgs, t2w_sin_wgs, & 1027 ft_max_error, cosft_duality_error, ierr) 1028 ABI_CHECK(ierr == 0, "Error in gx_minimax_grid") 1029 !if (ierr /= 0) then 1030 ! call gx_get_error_message(msg) 1031 ! ABI_ERROR(msg) 1032 !end if 1033 1034 call gaps%free() 1035 do iomega=1,Ep%nomegaei 1036 Ep%omega(Ep%nomegaer + iomega) = CMPLX(zero, iw_mesh(iomega), kind=dpc) 1037 !write(std_out, *)"iomega", Ep%omega(Ep%nomegaer + iomega) 1038 end do 1039 end if 1040 1041 else if (Ep%contour_deformation .and. Dtset%cd_customnimfrqs /= 0) then 1042 Ep%omega(Ep%nomegaer+1)=CMPLX(zero,Dtset%cd_imfrqs(1)) 1043 do iomega=2,Ep%nomegaei 1044 if (Dtset%cd_imfrqs(iomega) <= Dtset%cd_imfrqs(iomega-1)) then 1045 ABI_ERROR(' Specified imaginary frequencies need to be strictly increasing!') 1046 end if 1047 Ep%omega(Ep%nomegaer+iomega) = CMPLX(zero, Dtset%cd_imfrqs(iomega)) 1048 end do 1049 1050 else if (Ep%contour_deformation .and. Dtset%gw_frqim_inzgrid /= 0) then 1051 e0 = Dtset%ppmfrq; if (e0 < 0.1d-4) e0 = omegaplasma 1052 domegareal=one/(Ep%nomegaei+1) 1053 do iomega=1,Ep%nomegaei 1054 factor = iomega*domegareal 1055 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0*factor/(one-factor),kind=dpc) 1056 end do 1057 1058 else if (Ep%contour_deformation.and. Ep%nomegaei /= 0) then 1059 e0=Dtset%ppmfrq; if (e0<0.1d-4) e0=omegaplasma 1060 do iomega=1,Ep%nomegaei 1061 Ep%omega(Ep%nomegaer+iomega)=CMPLX(zero,e0/(Dtset%freqim_alpha-two)& 1062 * (EXP(two/(Ep%nomegaei+1)*LOG(Dtset%freqim_alpha-one)*iomega)-one),kind=dpc) 1063 end do 1064 end if 1065 1066 if (Dtset%cd_full_grid/=0) then 1067 ! Full grid will be calculated 1068 ! Grid values are added after the last imaginary freq. 1069 do ios=1,Ep%nomegaei 1070 do iomega=2,Ep%nomegaer 1071 Ep%omega(Ep%nomegaer+Ep%nomegaei+(ios-1)*(Ep%nomegaer-1)+(iomega-1)) = & 1072 CMPLX(REAL(Ep%omega(iomega)),AIMAG(Ep%omega(Ep%nomegaer+ios))) 1073 end do 1074 end do 1075 end if 1076 1077 ! Report frequency mesh for chi0. 1078 write(msg,'(2a)')ch10,' calculating chi0 at frequencies [eV] :' 1079 call wrtout([ab_out, std_out], msg) 1080 do iomega=1,Ep%nomega 1081 write(msg,'(i3,2es16.6)')iomega,Ep%omega(iomega)*Ha_eV 1082 call wrtout([ab_out, std_out], msg) 1083 end do 1084 1085 ! Allocate chi0, wings and array for chi0_sumrule check. 1086 ABI_MALLOC(chi0_sumrule, (Ep%npwe)) 1087 1088 write(msg,'(a,f12.1,a)')' Memory required for chi0 matrix= ',two*gwpc*Ep%npwe**2*Ep%nI*Ep%nJ*Ep%nomega*b2Mb," [Mb]." 1089 call wrtout(std_out, msg) 1090 ABI_MALLOC_OR_DIE(chi0, (Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr) 1091 ! 1092 !============================== END OF THE INITIALIZATION PART =========================== 1093 ! 1094 !====================================================================== 1095 !==== Loop over q-points. Calculate \epsilon^{-1} and save on disc ==== 1096 !====================================================================== 1097 call timab(321,2,tsec) ! screening(2) 1098 1099 iqcalc = 0 1100 if(Dtset%plowan_compute >= 10) then 1101 call init_plowannier(Dtset%plowan_bandf,Dtset%plowan_bandi,Dtset%plowan_compute,& 1102 Dtset%plowan_iatom,Dtset%plowan_it,Dtset%plowan_lcalc,Dtset%plowan_natom,& 1103 Dtset%plowan_nbl,Dtset%plowan_nt,Dtset%plowan_projcalc,Dtset%acell_orig,& 1104 Dtset%kptns,Dtset%nimage,Dtset%nkpt,Dtset%nspinor,Dtset%nsppol,Dtset%wtk,Dtset%dmft_t2g,wanibz_in) 1105 call get_plowannier(wanibz_in,wanibz,Dtset) 1106 call fullbz_plowannier(Dtset,Kmesh,Cryst,Pawang,wanibz,wanbz) 1107 end if 1108 1109 do iqibz=1,Qmesh%nibz 1110 call timab(306,1,tsec) 1111 is_first_qcalc=(iqibz==1) 1112 1113 ! Selective q-point calculation. 1114 found=.FALSE.; label=iqibz 1115 if (Ep%nqcalc /= Ep%nqibz) then 1116 do ii=1,Ep%nqcalc 1117 qtmp(:)=Qmesh%ibz(:,iqibz)-Ep%qcalc(:,ii) 1118 found=(normv(qtmp,gmet,'G')<GW_TOLQ) 1119 if (found) then 1120 label=ii; EXIT !ii 1121 end if 1122 end do 1123 if (.not.found) CYCLE !iqibz 1124 qtmp(:)=Ep%qcalc(:,1)-Qmesh%ibz(:,iqibz) 1125 is_first_qcalc=(normv(qtmp,gmet,'G')<GW_TOLQ) 1126 end if 1127 iqcalc = iqcalc + 1 1128 1129 bar=REPEAT('-',80) 1130 write(msg,'(4a,1x,a,i2,a,f9.6,2(",",f9.6),3a)')ch10,ch10,bar,ch10,& 1131 ' q-point number ',label,' q = (',(Qmesh%ibz(ii,iqibz),ii=1,3),') [r.l.u.]',ch10,bar 1132 call wrtout([ab_out, std_out], msg) 1133 is_qeq0 = 0; if (normv(Qmesh%ibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1 1134 1135 call timab(306,2,tsec) 1136 1137 if (is_qeq0 == 1) then 1138 ! Special treatment of the long wavelength limit. 1139 call timab(307,1,tsec) 1140 1141 ABI_MALLOC(chi0_head, (3,3,Ep%nomega)) 1142 ABI_MALLOC(chi0_lwing, (Ep%npwe*Ep%nI, Ep%nomega,3)) 1143 ABI_MALLOC(chi0_uwing, (Ep%npwe*Ep%nJ, Ep%nomega,3)) 1144 1145 call cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,qp_ebands,ks_ebands,Gsph_epsG0,& 1146 Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,nfftgw,& 1147 ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf,wanbz) 1148 1149 chihw = chi_new(ep%npwe, ep%nomega) 1150 chihw%head = chi0_head 1151 chihw%lwing = chi0_lwing 1152 chihw%uwing = chi0_uwing 1153 1154 ! Add the intraband term if required and metallic occupation scheme is used. 1155 add_chi0_intraband=.FALSE. !add_chi0_intraband=.TRUE. 1156 if (add_chi0_intraband .and. ebands_has_metal_scheme(qp_ebands)) then 1157 1158 ABI_MALLOC_OR_DIE(chi0intra,(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega), ierr) 1159 1160 ABI_MALLOC(chi0intra_lwing,(Ep%npwe*Ep%nI,Ep%nomega,3)) 1161 ABI_MALLOC(chi0intra_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,3)) 1162 ABI_MALLOC(chi0intra_head,(3,3,Ep%nomega)) 1163 1164 call chi0q0_intraband(Wfd,Cryst,Ep,Psps,qp_ebands,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,Dtset%usepawu,& 1165 ngfft_gw,chi0intra,chi0intra_head,chi0intra_lwing,chi0intra_uwing) 1166 1167 call wrtout(std_out,"Head of chi0 and chi0_intra") 1168 do iomega=1,Ep%nomega 1169 write(std_out,*)Ep%omega(iomega)*Ha_eV,REAL(chi0(1,1,iomega)),REAL(chi0intra(1,1,iomega)) 1170 write(std_out,*)Ep%omega(iomega)*Ha_eV,AIMAG(chi0(1,1,iomega)),AIMAG(chi0intra(1,1,iomega)) 1171 end do 1172 1173 chi0 = chi0 + chi0intra 1174 chi0_head = chi0_head + chi0intra_head 1175 chi0_lwing = chi0_lwing + chi0intra_lwing 1176 chi0_uwing = chi0_uwing + chi0intra_uwing 1177 1178 ABI_FREE(chi0intra) 1179 ABI_FREE(chi0intra_lwing) 1180 ABI_FREE(chi0intra_uwing) 1181 ABI_FREE(chi0intra_head) 1182 end if 1183 1184 if (.False.) then 1185 lwl_fname = strcat(dtfil%filnam_ds(4), "_LWL") 1186 call lwl_write(lwl_fname,cryst,vcp,ep%npwe,ep%nomega,gsph_epsg0%gvec,chi0,chi0_head,chi0_lwing,chi0_uwing,comm) 1187 end if 1188 1189 call timab(307,2,tsec) 1190 1191 else 1192 ! Calculate cchi0 for q/=0. 1193 call timab(308,1,tsec) 1194 call cchi0(use_tr,Dtset,Cryst,Qmesh%ibz(:,iqibz),Ep,Psps,Kmesh,qp_ebands,Gsph_epsG0,& 1195 Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftgw,ngfftf,nfftf_tot,chi0,ktabr,ktabrf,& 1196 Ltg_q(iqibz),chi0_sumrule,Wfd,Wfdf,wanbz) 1197 call timab(308,2,tsec) 1198 end if 1199 1200 ! Print chi0(q,G,Gp,omega), then calculate epsilon and epsilon^-1 for this q-point. 1201 ! Only master works but this part could be parallelized over frequencies. 1202 call timab(309,1,tsec) 1203 1204 do iomega=1,MIN(Ep%nomega, NOMEGA_PRINTED) 1205 write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 1206 call wrtout([ab_out, std_out], msg) 1207 write(msg,'(1x,a,i3,a,i4,a)')' chi0(q =',iqibz, ', omega =',iomega,', G,G'')' 1208 if (Ep%nqcalc /= Ep%nqibz) write(msg,'(a,i3,a,i4,a)')' chi0(q=',iqcalc,', omega=',iomega,', G,G'')' 1209 call wrtout(std_out, msg) 1210 ! arr99 is needed to avoid the update of all the tests. Now chi0 is divided by ucvol inside (cchi0|cchi0q0). 1211 ! TODO should be removed but GW tests have to be updated. 1212 ii = MIN(9, Ep%npwe) 1213 ABI_MALLOC(arr_99,(ii, ii)) 1214 arr_99 = chi0(1:ii,1:ii,iomega) * ucvol 1215 call print_arr(arr_99, max_r=2, unit=ab_out) 1216 call print_arr(arr_99, unit=std_out) 1217 ABI_FREE(arr_99) 1218 !call print_arr(chi0(:,:,iomega),max_r=2,unit=ab_out) 1219 !call print_arr(chi0(:,:,iomega),unit=std_out) 1220 end do 1221 1222 if (Ep%nomega > NOMEGA_PRINTED) then 1223 write(msg,'(a,i3,a)')' No. of calculated frequencies > ',NOMEGA_PRINTED,', stop printing ' 1224 call wrtout([ab_out, std_out], msg) 1225 end if 1226 1227 ! Write chi0 to _SUSC file 1228 ! Master creates and write the header if this is the first q-point calculated. 1229 if (Dtset%prtsuscep > 0 .and. my_rank == master) then 1230 title(1)='CHI0 file: chi0' 1231 title(2)=' ' 1232 if (is_qeq0 == 1) then 1233 string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string) 1234 title(1)=title(1)(1:21)//', calculated using inclvkb = '//string 1235 end if 1236 1237 ! Open file and write header for polarizability files. 1238 if (is_first_qcalc) then 1239 ikxc=0; test_type=0; tordering=1 1240 hchi0 = hscr_new("polarizability",dtset,ep,hdr_local,ikxc,test_type,tordering,title,Ep%npwe,Gsph_epsG0%gvec) 1241 1242 if (dtset%iomode == IO_MODE_ETSF) then 1243 NCF_CHECK(nctk_open_create(unt_susc, nctk_ncify(dtfil%fnameabo_sus), xmpi_comm_self)) 1244 NCF_CHECK(cryst%ncwrite(unt_susc)) 1245 NCF_CHECK(ebands_ncwrite(qp_ebands, unt_susc)) 1246 else 1247 unt_susc=Dtfil%unchi0 1248 if (open_file(dtfil%fnameabo_sus,msg,unit=unt_susc,status='unknown',form='unformatted') /= 0) then 1249 ABI_ERROR(msg) 1250 end if 1251 end if 1252 1253 fform_chi0 = hchi0%fform 1254 call hscr_io(hchi0,fform_chi0,2,unt_susc,xmpi_comm_self,0,Dtset%iomode) 1255 call Hchi0%free() 1256 end if 1257 1258 call write_screening("polarizability", unt_susc, Dtset%iomode, Ep%npwe, Ep%nomega, iqcalc, chi0) 1259 1260 if (dtset%iomode == IO_MODE_ETSF .and. is_qeq0 == 1 .and. Ep%nI == 1 .and. Ep%nJ == 1) then 1261 ! Write head and wings to file. See cchi0 for the equations needed to build chi0(q) for q--> 0. 1262 wing_shape = "two, three, number_of_coefficients_dielectric_function, number_of_frequencies_dielectric_function" 1263 ncerr = nctk_def_arrays(unt_susc, [ & 1264 nctkarr_t("sus_head", "dp", "two, three, three, number_of_frequencies_dielectric_function"), & 1265 nctkarr_t("sus_upper_wing", "dp", wing_shape), & 1266 nctkarr_t("sus_lower_wing", "dp", wing_shape) & 1267 ], defmode=.True.) 1268 NCF_CHECK(ncerr) 1269 1270 NCF_CHECK(nctk_set_datamode(unt_susc)) 1271 NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_head"), c2r(chi0_head))) 1272 1273 ABI_MALLOC(rwork_wing, (2, 3, Ep%npwe * Ep%nI, Ep%nomega)) 1274 do iomega=1,Ep%nomega 1275 do ig=1, Ep%npwe * Ep%nI 1276 rwork_wing(1,:,ig,iomega) = real(chi0_lwing(ig,iomega,:)) 1277 rwork_wing(2,:,ig,iomega) = aimag(chi0_lwing(ig,iomega,:)) 1278 end do 1279 end do 1280 NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_lower_wing"), rwork_wing)) 1281 1282 do iomega=1,Ep%nomega 1283 do ig=1, Ep%npwe * Ep%nI 1284 rwork_wing(1,:,ig,iomega) = real(chi0_uwing(ig,iomega,:)) 1285 rwork_wing(2,:,ig,iomega) = aimag(chi0_uwing(ig,iomega,:)) 1286 end do 1287 end do 1288 NCF_CHECK(nf90_put_var(unt_susc, nctk_idname(unt_susc, "sus_upper_wing"), rwork_wing)) 1289 ABI_FREE(rwork_wing) 1290 end if 1291 1292 end if ! is_first_qcalc 1293 1294 ! Calculate the Galitskii-Migdal and RPA functionals for the correlation energy if the polarizability on a 1295 ! Gauss-Legendre mesh along imaginary axis is available 1296 if (Ep%analytic_continuation .and. Dtset%gwrpacorr>0 ) then 1297 if (is_first_qcalc) then 1298 ABI_MALLOC(ec_rpa,(Dtset%gwrpacorr)) 1299 ec_rpa(:)=zero 1300 ec_gm=zero 1301 end if 1302 call calc_rpa_functional(Dtset%gwrpacorr,Dtset%gwgmcorr,label,iqibz,Ep,Vcp,Qmesh,Dtfil,gmet,chi0,comm,ec_rpa,ec_gm) 1303 if (label==Ep%nqcalc) then 1304 ABI_FREE(ec_rpa) 1305 end if 1306 end if 1307 1308 ! ========================================================== 1309 ! === Calculate RPA \tilde\epsilon^{-1} overwriting chi0 === 1310 ! ========================================================== 1311 approx_type=0 ! RPA 1312 option_test=0 ! TESTPARTICLE 1313 dim_wing=0; if (is_qeq0==1) dim_wing=3 1314 1315 if (dim_wing==0) then 1316 dim_wing=1 1317 if (.not.allocated(chi0_lwing)) then 1318 ABI_MALLOC(chi0_lwing,(Ep%npwe*Ep%nI,Ep%nomega,dim_wing)) 1319 end if 1320 if (.not.allocated(chi0_uwing)) then 1321 ABI_MALLOC(chi0_uwing,(Ep%npwe*Ep%nJ,Ep%nomega,dim_wing)) 1322 end if 1323 if (.not.allocated(chi0_head )) then 1324 ABI_MALLOC(chi0_head,(dim_wing,dim_wing,Ep%nomega)) 1325 end if 1326 dim_wing=0 1327 end if 1328 1329 #if 0 1330 ! Using the random q for the optical limit is one of the reasons 1331 ! why sigma breaks the initial energy degeneracies. 1332 chi0_lwing=czero 1333 chi0_uwing=czero 1334 chi0_head=czero 1335 #endif 1336 1337 ! Setup flags for the computation of em1 1338 ! If the vertex is being included for the spectrum, calculate the kernel now and pass it on 1339 if (dtset%gwgamma>0) rhor_kernel = rhor 1340 1341 select case (dtset%gwgamma) 1342 case (0) 1343 approx_type=0; option_test=0; dim_kxcg=0 1344 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1345 1346 case (1, 2) 1347 ! ALDA TDDFT kernel vertex 1348 ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1349 ABI_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!') 1350 ikxc=7; approx_type=1; dim_kxcg=1 1351 if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1352 if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1353 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1354 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,& 1355 Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self) 1356 1357 case (3, 4) 1358 ! ADA non-local kernel vertex 1359 ABI_CHECK(Wfd%usepaw==0,"ADA vertex + PAW not available") 1360 ABI_CHECK(Wfd%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases") 1361 ABI_WARNING('EXPERIMENTAL: Kernel is being added to screening, the SCR file will be non-standard!!') 1362 ikxc=7; approx_type=2 1363 if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1364 if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1365 ABI_MALLOC(fxc_ADA,(Ep%npwe,Ep%npwe,Ep%nqibz)) 1366 ! Use userrd to set kappa 1367 if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp 1368 ! Set correct value of kappa (should be scaled with alpha*r_s where) 1369 ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3) 1370 rhoav = (omegaplasma*omegaplasma)/four_pi 1371 r_s = (three/(four_pi*rhoav))**third 1372 alpha = (four*ninth*piinv)**third 1373 Dtset%userrd = Dtset%userrd*alpha*r_s 1374 1375 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf,Wfd%nspden,rhor_kernel,Ep%npwe,Ep%nqibz,Ep%qibz,& 1376 fxc_ADA,Gsph_epsG0%gvec,xmpi_comm_self,kappa_init=Dtset%userrd) 1377 1378 dim_kxcg=0 1379 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1380 1381 case (-3, -4, -5, -6, -7, -8) 1382 ! Bootstrap kernel and variants 1383 ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available") 1384 if (Dtset%gwgamma>-5) then 1385 ABI_WARNING('EXPERIMENTAL: Bootstrap kernel is being added to screening') 1386 approx_type=4 1387 else if (Dtset%gwgamma>-7) then 1388 ABI_WARNING('EXPERIMENTAL: Bootstrap kernel (head-only) is being added to screening') 1389 approx_type=5 1390 else 1391 ABI_WARNING('EXPERIMENTAL: RPA Bootstrap kernel is being added to screening') 1392 approx_type=6 1393 end if 1394 dim_kxcg=0 1395 option_test=MOD(Dtset%gwgamma,2) 1396 ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma 1397 ! 0 -> TESTPARTICLE, vertex in chi0 only 1398 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1399 1400 case (-11) 1401 ! LR+ALDA hybrid vertex kernel 1402 ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available") 1403 ikxc=7; dim_kxcg=1 1404 ABI_WARNING('EXPERIMENTAL: LR+ALDA hybrid kernel is being added to screening') 1405 approx_type=7 1406 option_test=1 ! TESTELECTRON 1407 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1408 rhor_kernel = rhor 1409 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,rhor_kernel,& 1410 Ep%npwe,dim_kxcg,kxcg,Gsph_epsG0%gvec,xmpi_comm_self) 1411 rhoav = (omegaplasma*omegaplasma)/four_pi 1412 1413 case default 1414 ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma))) 1415 end select 1416 1417 if (approx_type<2) then !ALDA 1418 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1419 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1420 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm) 1421 else if (approx_type<3) then !ADA 1422 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1423 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1424 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz)) 1425 else if (approx_type<7) then !Bootstrap 1426 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1427 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1428 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm) 1429 else if (approx_type<8) then !LR+ALDA 1430 call make_epsm1_driver(iqibz,dim_wing,Ep%npwe,Ep%nI,Ep%nJ,Ep%nomega,Ep%omega,& 1431 approx_type,option_test,Vcp,nfftf_tot,ngfftf,dim_kxcg,kxcg,Gsph_epsG0%gvec,& 1432 chi0_head,chi0_lwing,chi0_uwing,chi0,spectra,comm,rhor=rhoav) 1433 else 1434 ABI_ERROR(sjoin("Wrong approx_type:", itoa(approx_type))) 1435 end if 1436 1437 ABI_FREE(chi0_lwing) 1438 ABI_FREE(chi0_uwing) 1439 ABI_FREE(chi0_head) 1440 1441 if (my_rank == master .and. is_qeq0==1) then 1442 call spectra%repr(msg) 1443 call wrtout([ab_out, std_out], msg) 1444 1445 if (Ep%nomegaer > 2) then 1446 call spectra%write(W_EELF ,Dtfil%fnameabo_eelf) 1447 call spectra%write(W_EM_LF ,Dtfil%fnameabo_em1_lf) 1448 call spectra%write(W_EM_NLF,Dtfil%fnameabo_em1_nlf) 1449 end if 1450 end if ! master and is_qeq0==1 1451 1452 if (is_qeq0==1) call chi_free(chihw) 1453 1454 call spectra%free() 1455 ABI_SFREE(kxcg) 1456 ABI_SFREE(fxc_ADA) 1457 ! 1458 ! Output the sum rule evaluation. 1459 ! Vcp%vc_sqrt(:,iqibz) Contains vc^{1/2}(q,G), complex-valued due to a possible cutoff 1460 epsm1 => chi0 1461 call output_chi0sumrule((is_qeq0==1),iqibz,Ep%npwe,omegaplasma,chi0_sumrule,epsm1(:,:,1),Vcp%vc_sqrt(:,iqibz)) 1462 1463 ! If input variable npvel is larger than 0, trigger the Random Stopping Power calculation 1464 ! Only the masternode is used 1465 if (my_rank==master .and. Dtset%npvel>0) then 1466 if (is_first_qcalc) then 1467 ABI_MALLOC(rspower,(Dtset%npvel)) 1468 rspower(:)=zero 1469 end if 1470 call random_stopping_power(iqibz,Dtset%npvel,Dtset%pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower) 1471 if (label==Ep%nqcalc) then 1472 ABI_FREE(rspower) 1473 end if 1474 end if 1475 1476 ! Write heads and wings to main output file. 1477 if (is_qeq0 == 1) then 1478 write(msg,'(1x,2a)')' Heads and wings of the symmetrical epsilon^-1(G,G'') ',ch10 1479 call wrtout(ab_out,msg) 1480 do iomega=1,Ep%nomega 1481 write(msg,'(2x,a,i4,a,2f9.4,a)')' Upper and lower wings at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 1482 call wrtout(ab_out, msg) 1483 call print_arr(epsm1(1,:,iomega),max_r=9,unit=ab_out) 1484 call print_arr(epsm1(:,1,iomega),max_r=9,unit=ab_out) 1485 call wrtout(ab_out, ch10) 1486 end do 1487 end if 1488 1489 call timab(309,2,tsec) 1490 call timab(310,1,tsec) ! wrscr 1491 1492 if (my_rank==master) then 1493 ! === Write the symmetrical epsilon^-1 on file === 1494 title(1)='SCR file: epsilon^-1' 1495 if (is_qeq0==1) then 1496 string='0'; if (Dtset%usepaw==0.and.Ep%inclvkb/=0) call int2char10(Ep%inclvkb,string) 1497 title(1)=title(1)(1:21)//', calculated using inclvkb = '//string 1498 end if 1499 title(2)='TESTPARTICLE' 1500 ctype='RPA' 1501 title(2)(14:17)=ctype !this has to be modified 1502 1503 if (is_first_qcalc) then 1504 ! === Open file and write the header for the SCR file === 1505 ! * Here we write the RPA approximation for \tilde\epsilon^{-1} 1506 ikxc=0; test_type=0; tordering=1 1507 hem1 = hscr_new("inverse_dielectric_function",dtset,ep,hdr_local,ikxc,test_type,tordering,title,& 1508 & Ep%npwe,Gsph_epsG0%gvec) 1509 fform_em1 = hem1%fform 1510 if (dtset%iomode == IO_MODE_ETSF) then 1511 NCF_CHECK(nctk_open_create(unt_em1, nctk_ncify(dtfil%fnameabo_scr), xmpi_comm_self)) 1512 NCF_CHECK(cryst%ncwrite(unt_em1)) 1513 NCF_CHECK(ebands_ncwrite(qp_ebands, unt_em1)) 1514 else 1515 unt_em1=Dtfil%unscr 1516 if (open_file(dtfil%fnameabo_scr,msg,unit=unt_em1,status='unknown',form='unformatted') /= 0) then 1517 ABI_ERROR(msg) 1518 end if 1519 end if 1520 call hscr_io(hem1,fform_em1,2,unt_em1,xmpi_comm_self,0,Dtset%iomode) 1521 call Hem1%free() 1522 end if 1523 1524 call write_screening("inverse_dielectric_function",unt_em1,Dtset%iomode,Ep%npwe,Ep%nomega,iqcalc,epsm1) 1525 end if ! my_rank==master 1526 1527 call timab(310,2,tsec) 1528 end do ! Loop over q-points 1529 1530 if (Dtset%plowan_compute >= 10) call destroy_plowannier(wanbz) 1531 1532 ! Close Files. 1533 if (my_rank == master) then 1534 if (dtset%iomode == IO_MODE_ETSF) then 1535 NCF_CHECK(nf90_close(unt_em1)) 1536 if (dtset%prtsuscep > 0) then 1537 NCF_CHECK(nf90_close(unt_susc)) 1538 end if 1539 else 1540 close(unt_em1) 1541 if (dtset%prtsuscep > 0) close(unt_susc) 1542 end if 1543 end if 1544 ! 1545 !===================== 1546 !==== Free memory ==== 1547 !===================== 1548 ABI_FREE(chi0_sumrule) 1549 ABI_FREE(chi0) 1550 ABI_SFREE(rhor_kernel) 1551 1552 ABI_FREE(rhor) 1553 ABI_FREE(rhog) 1554 ABI_FREE(ks_vbik) 1555 ABI_FREE(qp_vbik) 1556 ABI_FREE(ktabr) 1557 ABI_FREE(taur) 1558 ABI_FREE(ks_vhartr) 1559 ABI_FREE(ks_vtrial) 1560 ABI_FREE(vpsp) 1561 ABI_FREE(ks_vxc) 1562 ABI_FREE(ph1d) 1563 ABI_FREE(ph1df) 1564 ABI_FREE(nhatgr) 1565 ABI_FREE(nhat) 1566 call pawfgr_destroy(Pawfgr) 1567 1568 if (Dtset%usepaw==1) then ! Optional deallocation for PAW. 1569 call pawrhoij_free(Pawrhoij) 1570 call pawfgrtab_free(Pawfgrtab) 1571 call paw_ij_free(Paw_ij) 1572 call paw_an_free(Paw_an) 1573 call pawpwff_free(Paw_pwff) 1574 if (Dtset%pawcross==1) then 1575 call paw_pwaves_lmn_free(Paw_onsite) 1576 call wfdf%free() 1577 end if 1578 end if 1579 1580 ABI_FREE(Pawfgrtab) 1581 ABI_FREE(Paw_pwff) 1582 ABI_FREE(Pawrhoij) 1583 ABI_FREE(Paw_ij) 1584 ABI_FREE(Paw_an) 1585 ABI_FREE(ktabrf) 1586 ABI_FREE(Paw_onsite) 1587 1588 call wfd%free() 1589 call Kmesh%free() 1590 call Qmesh%free() 1591 call cryst%free() 1592 call Gsph_epsG0%free() 1593 call Gsph_wfn%free() 1594 call Vcp%free() 1595 call Ep%free() 1596 call Hdr_wfk%free() 1597 call Hdr_local%free() 1598 call ebands_free(ks_ebands) 1599 call ebands_free(qp_ebands) 1600 call destroy_mpi_enreg(MPI_enreg_seq) 1601 call littlegroup_free(ltg_q) 1602 ABI_FREE(Ltg_q) 1603 1604 call timab(301,2,tsec) 1605 1606 DBG_EXIT("COLL") 1607 1608 end subroutine screening
m_screening_driver/setup_screening [ Functions ]
[ Top ] [ m_screening_driver ] [ Functions ]
NAME
setup_screening
FUNCTION
Initialize the Ep% data type containing the parameters used during the screening calculation. as well as basic objects describing the BZ sampling .... TODO list to be completed
INPUTS
wfk_fname=Name of the input WFK file. acell(3)=length scales of primitive translations (Bohr). rprim(3,3)=dimensionless real space primitive translations. dtfil <type(datafiles_type)>=variables related to files
OUTPUT
ngfft_gw(18)=Contain all needed information about the 3D FFT for the oscillator strengths. See ~abinit/doc/variables/vargs.htm#ngfft Ltg_q(:)<littlegroup_t>,= Ep<em1params_t>=Parameters for the screening calculation. Most part of it is Initialized and checked. Hdr_wfk type(Hdr_type)=Header of the KSS file. Cryst<crystal_t>=Definition of the unit cell and its symmetries. Kmesh<kmesh_t>=Structure defining the k-point sampling (wavefunctions). Qmesh<kmesh_t>=Structure defining the q-point sampling (screening) Gsph_wfn<gsphere_t>=Structure defining the G-sphere for the wavefunctions (not k-dependent). Gsph_epsG0<gsphere_t>=The G-sphere for the screening, enlarged to take into account for umklapps. Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file Vcp <type vcoul_t> datatype gathering information on the coulombian cutoff technique comm=MPI communicator.
SIDE EFFECTS
Dtset<Dataset_type>=All input variables for this dataset. %ecutwfn, %npwwfn, %ecuteps, %npweps might be redefined in setshells in order to close the shell.
SOURCE
1649 subroutine setup_screening(codvsn,acell,rprim,wfk_fname,Dtset,Psps,Pawtab,& 1650 ngfft_gw,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,ks_ebands,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm) 1651 1652 !Arguments ------------------------------------ 1653 !scalars 1654 integer,intent(in) :: comm 1655 character(len=8),intent(in) :: codvsn 1656 character(len=fnlen),intent(in) :: wfk_fname 1657 type(Dataset_type),intent(inout) :: Dtset !INOUT is due to setshells 1658 type(Pseudopotential_type),intent(in) :: Psps 1659 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 1660 type(em1params_t),intent(out) :: Ep 1661 type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out 1662 type(ebands_t),intent(out) :: ks_ebands 1663 type(kmesh_t),intent(out) :: Kmesh,Qmesh 1664 type(crystal_t),intent(out) :: Cryst 1665 type(gsphere_t),intent(out) :: Gsph_epsG0,Gsph_wfn 1666 type(vcoul_t),intent(out) :: Vcp 1667 !arrays 1668 integer,intent(out) :: ngfft_gw(18) 1669 real(dp),intent(in) :: acell(3),rprim(3,3) 1670 type(littlegroup_t),pointer :: Ltg_q(:) 1671 1672 !Local variables------------------------------- 1673 !scalars 1674 integer,parameter :: NOMEGAGAUSS=30,NOMEGAREAL=201,pertcase0=0,master=0 1675 integer :: bantot,ib,ibtot,ikibz,iq,iqp,isppol,ig,ng,ierr 1676 integer :: jj,mod10,mband,ng_kss,iqbz,isym,iq_ibz,itim 1677 integer :: timrev,use_umklp !,ncerr 1678 integer :: npwepG0,nshepspG0,method,enforce_sym,nfftgw_tot !,spin,band,ik_ibz, 1679 integer :: istart,iend,test_npwkss,my_rank,nprocs !ii 1680 real(dp),parameter :: OMEGAERMAX=100.0/Ha_eV 1681 real(dp) :: ecutepspG0,ucvol,domegareal 1682 logical :: remove_inv,ltest,found,is_static,has_q0 1683 character(len=500) :: msg 1684 type(wvl_internal_type) :: wvl 1685 !arrays 1686 integer :: ng0sh_opt(3) 1687 integer,allocatable :: npwarr(:) 1688 integer,pointer :: gvec_kss(:,:) 1689 integer,pointer :: test_gvec_kss(:,:) 1690 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),qtmp(3),sq(3),qbz(3) 1691 real(dp),pointer :: energies_p(:,:,:) 1692 real(dp),allocatable :: doccde(:),eigen(:),occfact(:) 1693 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 1694 1695 ! ************************************************************************* 1696 1697 DBG_ENTER('COLL') 1698 1699 ! Check for calculations that are not implemented 1700 ltest = ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol) == Dtset%nband(1)) 1701 ABI_CHECK(ltest, 'dtset%nband(:) must be constant in the GW code.') 1702 1703 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1704 1705 call mkrdim(acell,rprim,rprimd) 1706 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1707 1708 ! Set up basic parameters of the calculation 1709 Ep%gwcalctyp =Dtset%gwcalctyp 1710 Ep%plasmon_pole_model =.TRUE. 1711 Ep%analytic_continuation=.FALSE. 1712 Ep%contour_deformation =.FALSE. 1713 1714 mod10=MOD(Ep%gwcalctyp,10) 1715 if (mod10/=0.and.mod10/=8) Ep%plasmon_pole_model =.FALSE. 1716 if (mod10==1) Ep%analytic_continuation=.TRUE. 1717 if (mod10==2.or.mod10==9) Ep%contour_deformation =.TRUE. 1718 is_static=(mod10==5.or.mod10==6.or.mod10==7) 1719 1720 Ep%nbnds =Dtset%nband(1) 1721 Ep%symchi =Dtset%symchi 1722 Ep%inclvkb=Dtset%inclvkb; if (Dtset%usepaw/=0) Ep%inclvkb=0 1723 Ep%zcut =Dtset%zcut 1724 1725 write(msg,'(2a,i4,2a,f10.6,a)')ch10,& 1726 ' GW calculation type = ',Ep%gwcalctyp,ch10,& 1727 ' zcut to avoid poles in chi0 [eV] = ',Ep%zcut*Ha_eV,ch10 1728 call wrtout(std_out, msg) 1729 1730 Ep%awtr =Dtset%awtr 1731 Ep%npwe =Dtset%npweps 1732 Ep%npwwfn=Dtset%npwwfn 1733 Ep%npwvec=MAX(Ep%npwe,Ep%npwwfn) 1734 1735 timrev = 2 ! This information is not reported in the header 1736 ! 1 --> do not use time-reversal symmetry 1737 ! 2 --> take advantage of time-reversal symmetry 1738 if (any(dtset%kptopt == [3, 4])) timrev = 1 1739 1740 if (timrev==1.and.Dtset%awtr/=0) then 1741 ABI_ERROR("awtr/=0 cannot be used when time-reversal symmetry doesn't hold") 1742 end if 1743 1744 ! Read parameters from WFK and verifify them. 1745 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 1746 mband = MAXVAL(Hdr_wfk%nband) 1747 call hdr_wfk%vs_dtset(dtset) 1748 remove_inv=.FALSE. 1749 1750 test_npwkss = 0 1751 call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,& 1752 gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr) 1753 ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss") 1754 1755 ABI_MALLOC(gvec_kss,(3,test_npwkss)) 1756 gvec_kss = test_gvec_kss 1757 ng_kss = test_npwkss 1758 1759 if (Ep%npwvec>ng_kss) then 1760 Ep%npwvec=ng_kss 1761 if (Ep%npwwfn> ng_kss) Ep%npwwfn=ng_kss 1762 if (Ep%npwe > ng_kss) Ep%npwe =ng_kss 1763 write(msg,'(3a,3(a,i6,a))')ch10,& 1764 ' Number of G-vectors found less then required. Calculation will proceed with ',ch10,& 1765 ' npwvec = ',Ep%npwvec,ch10,& 1766 ' npweps = ',Ep%npwe ,ch10,& 1767 ' npwwfn = ',Ep%npwwfn,ch10 1768 ABI_WARNING(msg) 1769 end if 1770 1771 ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2)) 1772 ierr = 0 1773 do ig=1,ng 1774 if (ANY(gvec_kss(:,ig)/=test_gvec_kss(:,ig))) then 1775 ierr=ierr+1 1776 write(std_out,*)" gvec_kss ",ig,"/",ng,gvec_kss(:,ig),test_gvec_kss(:,ig) 1777 end if 1778 end do 1779 ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss") 1780 ABI_FREE(test_gvec_kss) 1781 1782 ! Get important dimension from Hdr_wfk 1783 ! Check also the consistency btw Hdr_wfk and Dtset. 1784 Ep%nsppol=Hdr_wfk%nsppol 1785 Ep%nkibz =Hdr_wfk%nkpt 1786 1787 if (Ep%nbnds>mband) then 1788 Ep%nbnds=mband 1789 Dtset%nband(:)=mband 1790 write(msg,'(4a,i4,a)')ch10,& 1791 ' Number of bands found less then required. ',ch10,& 1792 ' Calculation will proceed with nbnds = ',mband,ch10 1793 ABI_WARNING(msg) 1794 end if 1795 1796 cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv) 1797 call cryst%print(mode_paral='COLL') 1798 1799 ! === Create basic data types for the calculation === 1800 ! * Kmesh defines the k-point sampling for the wavefunctions. 1801 ! * Qmesh defines the q-point sampling for chi0, all possible differences k1-k2 reduced to the IBZ. 1802 ! TODO Kmesh%bz should be [-half,half[ but this modification will be painful! 1803 1804 call Kmesh%init(Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.) 1805 !call Kmesh%init(Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.) 1806 ! Some required information are not filled up inside kmesh_init 1807 ! So doing it here, even though it is not clean 1808 Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:) 1809 Kmesh%nshift =Dtset%nshiftk 1810 ABI_MALLOC(Kmesh%shift,(3,Kmesh%nshift)) 1811 Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk) 1812 1813 call Kmesh%print("K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 1814 call Kmesh%print("K-mesh for the wavefunctions",ab_out, 0, "COLL") 1815 1816 ! === Find Q-mesh, and do setup for long wavelength limit === 1817 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 1818 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 1819 call find_qmesh(Qmesh,Cryst,Kmesh) 1820 1821 call Qmesh%print("Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL") 1822 call Qmesh%print("Q-mesh for the screening function",ab_out ,0 ,"COLL") 1823 1824 do iqbz=1,Qmesh%nbz 1825 call qmesh%get_BZ_item(iqbz,qbz,iq_ibz,isym,itim) 1826 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 1827 if (ANY(ABS(qbz-sq )>1.0d-4)) then 1828 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 1829 ' qpoint ',qbz,' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 1830 ' through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 1831 ' however a non zero umklapp G_o vector is required and this is not yet allowed' 1832 ABI_ERROR(msg) 1833 end if 1834 end do 1835 1836 ! This section is now performed in invars2 1837 ! Write the list of qpoints for the screening in netcdf format and exit. 1838 ! This file is used by abipy to generate multiple input files. 1839 ! if (Dtset%nqptdm == -1) then 1840 ! if (my_rank==master) then 1841 ! ncerr = nctk_write_ibz(strcat(dtfil%filnam_ds(4), "_qptdms.nc"), qmesh%ibz, qmesh%wt) 1842 ! NCF_CHECK(ncerr) 1843 ! end if 1844 ! ABI_ERROR_NODUMP("Aborting now") 1845 ! end if 1846 1847 if (Dtset%gw_nqlwl==0) then 1848 Ep%nqlwl=1 1849 ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl)) 1850 Ep%qlwl(:,1)=GW_Q0_DEFAULT ! Use default shift 0.000010, 0.000020, 0.000030 1851 else 1852 Ep%nqlwl=Dtset%gw_nqlwl 1853 ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl)) 1854 Ep%qlwl(:,:)=Dtset%gw_qlwl(:,1:Ep%nqlwl) 1855 ABI_CHECK(Ep%nqlwl==1,"nqlwl/=1 not coded") 1856 end if 1857 !write(std_out,*)" Using qlwl = ",Ep%qlwl 1858 1859 ! Find optimal value for G-sphere enlargment due to oscillator matrix elements 1860 call get_ng0sh(Kmesh%nbz,Kmesh%bz,Qmesh%nibz,Qmesh%ibz,Kmesh%nbz,Kmesh%bz,GW_TOLQ0,ng0sh_opt) 1861 call wrtout(std_out,sjoin(' Optimal value for ng0sh:',ltoa(ng0sh_opt)),"COLL") 1862 1863 Ep%mG0(:)=ng0sh_opt(:) !Ep%mG0(:) = [3, 3, 3] 1864 1865 ! === In case of symmetrization, find the little group of the q"s === 1866 ! * For the long-wavelength limit we consider a small but finite q. However the oscillators are 1867 ! evaluated setting q==0. Thus it is possible to take advantage of symmetries also when q --> 0. 1868 ! * Here we calculate the enlargement of the G-sphere, npwepG0, needed to account for umklapps. 1869 ! TODO Switch on use_umklp, write all this stuff to ab_out 1870 1871 Ep%npwepG0 = Ep%npwe 1872 ABI_MALLOC(Ltg_q, (Qmesh%nibz)) 1873 1874 do iq=1,Qmesh%nibz 1875 qtmp = Qmesh%ibz(:,iq); if (normv(qtmp,gmet,'G') < GW_TOLQ0) qtmp(:) = zero; use_umklp = 0 1876 call Ltg_q(iq)%init(qtmp, Kmesh%nbz, Kmesh%bz, Cryst, use_umklp, Ep%npwe, gvec=gvec_kss) 1877 end do 1878 1879 ecutepspG0 = Dtset%ecuteps 1880 ABI_CHECK(ecutepspG0 > zero, "ecuteps must be > 0") 1881 if (Ep%symchi/=0) then 1882 ecutepspG0=MAXVAL(Ltg_q(:)%max_kin_gmG0)+tol6; npwepG0=0; nshepspG0=0 1883 if (my_rank == master) write(std_out,*)" Due to umklapp processes : ecutepspg0= ",ecutepspG0 1884 call setshells(ecutepspG0,npwepG0,nshepspG0,Cryst%nsym,gmet,gprimd,Cryst%symrel,'eps_pG0',Cryst%ucvol) 1885 Ep%npwepG0=npwepG0 1886 end if 1887 1888 if (Ep%npwepG0>Ep%npwvec) then 1889 write(msg,'(3a,i5,a,i5)')& 1890 ' npwepG0 > npwvec, decrease npweps or increase npwwfn. ',ch10,& 1891 ' npwepG0 = ',Ep%npwepG0,' npwvec = ',Ep%npwvec 1892 ABI_ERROR(msg) 1893 end if 1894 1895 ! === Create structure describing the G-sphere used for chi0/espilon and Wfns === 1896 ! * The cutoff is >= ecuteps to allow for umklapp 1897 #if 0 1898 call Gsph_wfn%init(Cryst, 0, ecut=Dtset%ecutwfn) 1899 Dtset%npwwfn = Gsph_wfn%ng 1900 Ep%npwwfn = Gsph_wfn%ng 1901 ierr = 0 1902 do ig=1,MIN(Gsph_wfn%ng, ng_kss) 1903 if ( ANY(Gsph_wfn%gvec(:,ig) /= gvec_kss(:,ig)) ) then 1904 write(std_out,*)ig, Gsph_wfn%gvec(:,ig), gvec_kss(:,ig) 1905 end if 1906 end do 1907 ABI_CHECK(ierr==0,"Wrong gvec_wfn") 1908 #else 1909 call Gsph_wfn%init(Cryst, Ep%npwvec, gvec=gvec_kss) 1910 #endif 1911 1912 #if 0 1913 call Gsph_epsG0%init(Cryst, 0, ecut=ecutepspG0) 1914 Ep%npwepG0 = Gsph_epsG0%ng 1915 ierr = 0 1916 do ig=1,MIN(Gsph_epsG0%ng, ng_kss) 1917 if ( ANY(Gsph_epsG0%gvec(:,ig) /= gvec_kss(:,ig)) ) then 1918 write(std_out,*)ig, Gsph_epsG0%gvec(:,ig), gvec_kss(:,ig) 1919 end if 1920 end do 1921 ABI_CHECK(ierr==0,"Wrong gvec_epsG0") 1922 #else 1923 call Gsph_epsG0%init(Cryst, Ep%npwepG0, gvec=gvec_kss) 1924 #endif 1925 ! 1926 ! ======================================================================= 1927 ! ==== Setup of the FFT mesh used for the oscillator matrix elements ==== 1928 ! ======================================================================= 1929 ! * ngfft_gw(7:18) is the same as Dtset%ngfft(7:18), initialized before entering setup_screening. 1930 ! Here we just redefine ngfft_gw(1:6) according to the following options: 1931 ! 1932 ! method==0 ==> FFT grid read from __fft.in__ (only for debugging purpose) 1933 ! method==1 ==> normal FFT grid 1934 ! method==2 ==> slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 1935 ! method==3 ==> doubled FFT grid, to treat exactly the convolution defining the density, 1936 ! Useful in sigma if ppmodel=[2,3,4] since rho(G-Gp) or to calculate matrix elements of v_Hxc. 1937 ! 1938 ! enforce_sym==1 ==> enforce a direct space FFT mesh compatible with all symmetries operation 1939 ! enforce_sym==0 ==> Find the smallest FFT grid compatibile with the library, do not care about symmetries 1940 ! 1941 ngfft_gw(1:18)=Dtset%ngfft(1:18); method=2 1942 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 1943 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 1944 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 1945 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 1946 enforce_sym=MOD(Dtset%fftgw,10) 1947 1948 ! Use npwepG0 to account for umklapps. 1949 call setmesh(gmet,gvec_kss,ngfft_gw,Ep%npwvec,Ep%npwepG0,Ep%npwwfn,nfftgw_tot,method,Ep%mG0,Cryst,enforce_sym) 1950 !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Ep%mG0,enforce_sym,ngfft_gw,nfftgw_tot) 1951 1952 ABI_FREE(gvec_kss) 1953 1954 ! FIXME this wont work if nqptdm/=0 1955 call Vcp%init(Gsph_epsG0,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecuteps,Ep%npwe,Ep%nqlwl,& 1956 Ep%qlwl,comm) 1957 1958 #if 0 1959 ! Using the random q for the optical limit is one of the reasons 1960 ! why sigma breaks the initial energy degeneracies. 1961 Vcp%i_sz=zero 1962 Vcp%vc_sqrt(1,:)=czero 1963 Vcp%vcqlwl_sqrt(1,:)=czero 1964 #endif 1965 1966 ! Value of scissor energy 1967 Ep%mbpt_sciss=Dtset%mbpt_sciss 1968 1969 ! Define the frequency mesh for epsilon according to the method used. 1970 Ep%nomegaei=1 1971 Ep%nomegaer=1; if (is_static) Ep%nomegaer=0 1972 Ep%nomegaec=0 1973 Ep%omegaermax=zero 1974 1975 ! For ppmodels 2,3,4, only omega=0 is needed. 1976 if (Ep%plasmon_pole_model.and.Dtset%nfreqre==1.and.Dtset%nfreqim==0) then 1977 Ep%nomegaer=1; Ep%nomegaei=0 1978 write(msg,'(7a)')ch10,& 1979 ' The inverse dielectric matrix will be calculated on zero frequency only',ch10,& 1980 ' please note that the calculated epsilon^-1 cannot be used ',ch10,& 1981 ' to calculate QP corrections using plasmonpole model 1',ch10 1982 call wrtout([ab_out, std_out], msg) 1983 end if 1984 1985 ! Max number of omega along the imaginary axis 1986 if (Ep%analytic_continuation.or.Ep%contour_deformation) then 1987 Ep%nomegaei=Dtset%nfreqim 1988 if (Dtset%gw_frqim_inzgrid==1) then 1989 ABI_WARNING('iomega = z/1-z transfom grid will be used for imaginary frequency grid') 1990 end if 1991 if (Dtset%cd_customnimfrqs/=0) then 1992 ABI_WARNING('Custom imaginary grid specified. Assuming experienced user.') 1993 Ep%nomegaei=Dtset%cd_customnimfrqs 1994 end if 1995 if (Ep%nomegaei==-1) then 1996 Ep%nomegaei=NOMEGAGAUSS 1997 ABI_WARNING(sjoin('Number of imaginary frequencies set to default= ',itoa(NOMEGAGAUSS))) 1998 end if 1999 if (Ep%nomegaei==0) then 2000 ABI_WARNING(' nfreqim = 0! Assuming experienced user merging several frequency calculations.') 2001 end if 2002 end if 2003 2004 ! Range and total number of real frequencies. 2005 Ep%omegaermin = zero 2006 if (Ep%contour_deformation) then 2007 Ep%nomegaer=Dtset%nfreqre; Ep%omegaermin=Dtset%freqremin; Ep%omegaermax=Dtset%freqremax 2008 if (Dtset%gw_frqre_tangrid==1) then 2009 Ep%omegaermax=Dtset%cd_max_freq 2010 ABI_WARNING('Tangent transfom grid will be used for real frequency grid') 2011 end if 2012 if (Dtset%gw_frqre_tangrid==1) then 2013 ABI_WARNING('Tangent transfom grid will be used for real frequency grid') 2014 end if 2015 if (Ep%nomegaer==-1) then 2016 Ep%nomegaer=NOMEGAREAL 2017 ABI_WARNING(sjoin('Number of real frequencies set to default= ',itoa(NOMEGAREAL))) 2018 end if 2019 if (Ep%nomegaer==0) then 2020 ABI_WARNING('nfreqre = 0 ! Assuming experienced user.') 2021 end if 2022 if (ABS(Ep%omegaermin)<TOL16) then 2023 Ep%omegaermin=zero 2024 write(msg,'(a,f8.4)')' Min real frequency set to default [Ha] = ',Ep%omegaermin 2025 ABI_WARNING(msg) 2026 end if 2027 if (Ep%omegaermin>Ep%omegaermax) then 2028 ABI_ERROR('freqremin > freqremax !') 2029 end if 2030 if (Ep%omegaermax<TOL16) then 2031 Ep%omegaermax=OMEGAERMAX 2032 write(msg,'(a,f8.4)')' Max real frequency set to default [Ha] = ',OMEGAERMAX 2033 ABI_WARNING(msg) 2034 end if 2035 ! Check if a subset of the frequencies is to be used 2036 if (Dtset%cd_subset_freq(1)/=0) then 2037 istart = Dtset%cd_subset_freq(1) 2038 iend = Dtset%cd_subset_freq(2) 2039 if (istart>iend.or.istart<0.or.iend<0) then 2040 ABI_ERROR(' check indices of cd_subset_freq!') 2041 end if 2042 write(msg,'(2(a,i0))')' Using cd_subset_freq to only do freq. from ',istart,' to ',iend 2043 ABI_WARNING(msg) 2044 ! Reset the numbers 2045 if (Dtset%gw_frqre_tangrid/=1) then ! Normal equidistant grid 2046 Ep%nomegaer = iend-istart+1 2047 domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1) 2048 Ep%omegaermin = Ep%omegaermin+(istart-1)*domegareal 2049 Ep%omegaermax = Ep%omegaermin+(iend-1)*domegareal 2050 else 2051 Ep%nomegaer = iend-istart+1 2052 end if 2053 end if 2054 end if 2055 2056 ! Check full grid calculations 2057 if (Dtset%cd_full_grid/=0) then 2058 ABI_WARNING("FULL GRID IN COMPLEX PLANE CALCULATED. YOU MIGHT NOT BE ABLE TO USE SCREENING FILES!") 2059 if (Dtset%cd_subset_freq(1)/=0) then 2060 ABI_ERROR('cd_subset_freq cannot be used with cd_full_grid!') 2061 end if 2062 Ep%nomegaec = Ep%nomegaei*(Ep%nomegaer-1) 2063 end if 2064 2065 Ep%nomega=Ep%nomegaer+Ep%nomegaei+Ep%nomegaec ! Total number of frequencies. 2066 2067 ! ==== Setup of the spectral method ==== 2068 Ep%spmeth =Dtset%spmeth; Ep%nomegasf=Dtset%nomegasf; Ep%spsmear =Dtset%spbroad 2069 2070 if (Ep%spmeth/=0) then 2071 write(msg,'(2a,i3,2a,i8)')ch10,& 2072 & ' setup_screening: using spectral method: ',Ep%spmeth,ch10,& 2073 & ' Number of frequencies for imaginary part: ',Ep%nomegasf 2074 call wrtout(std_out, msg) 2075 if (Ep%spmeth==2) then 2076 write(msg,'(a,f8.5,a)')' Gaussian broadening = ',Ep%spsmear*Ha_eV,' [eV]' 2077 call wrtout(std_out, msg) 2078 end if 2079 end if 2080 2081 Ep%nI=1; Ep%nJ=1 2082 if (Dtset%nspinor==2) then 2083 !if (Dtset%usepaw==1.and.Dtset%pawspnorb>0) then 2084 ! Ep%nI=1; Ep%nJ=4 2085 !end if 2086 ! For spin-spin interaction 2087 ! Ep%nI=4; Ep%nJ=4 2088 ABI_CHECK(Ep%npwepG0 == Ep%npwe, "npwepG0 must be == npwe if spinor==2") 2089 !ABI_CHECK(Ep%symchi == 0, "symchi/=0 and nspinor=2 not available") 2090 end if 2091 2092 ! === Enable the calculations of chi0 on user-specified q-points === 2093 Ep%nqibz=Qmesh%nibz 2094 ABI_MALLOC(Ep%qibz,(3,Ep%nqibz)) 2095 Ep%qibz(:,:)=Qmesh%ibz(:,:) 2096 2097 Ep%nqcalc=Ep%nqibz 2098 if (Dtset%nqptdm>0) Ep%nqcalc=Dtset%nqptdm 2099 2100 ABI_MALLOC(Ep%qcalc,(3,Ep%nqcalc)) 2101 if (Ep%nqcalc/=Ep%nqibz) then 2102 write(msg,'(6a)')ch10,& 2103 ' Dielectric matrix will be calculated only for some ',ch10,& 2104 ' selected q points provided by the user through the input variables ',ch10,& 2105 ' nqptdm and qptdm' 2106 call wrtout([ab_out, std_out], msg) 2107 ltest= Ep%nqcalc <= Qmesh%nibz 2108 ABI_CHECK(ltest, 'nqptdm should not exceed the number of q points in the IBZ') 2109 Ep%qcalc(:,:)=Dtset%qptdm(:,1:Ep%nqcalc) 2110 ! Check whether the q-points provided are correct. 2111 do iq=1,Ep%nqcalc 2112 found=.FALSE. 2113 do iqp=1,Qmesh%nibz 2114 qtmp(:)=Ep%qcalc(:,iq)-Qmesh%ibz(:,iqp) 2115 found=(normv(qtmp,gmet,'G')<GW_TOLQ) 2116 if (found) EXIT 2117 end do 2118 ABI_CHECK(found, 'One or more points specified by Dtset%qptdm do not satisfy q=k1-k2') 2119 end do 2120 else 2121 Ep%qcalc(:,:)=Ep%qibz(:,:) 2122 end if 2123 2124 ! To write the SCR header correctly, with heads and wings, we have 2125 ! to make sure that q==0, if present, is the first q-point in the list. 2126 has_q0=.FALSE. 2127 do iq=1,Ep%nqcalc 2128 if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then 2129 has_q0=.TRUE.; EXIT 2130 end if 2131 end do 2132 2133 if (has_q0.and.normv(Ep%qcalc(:,1),gmet,'G')>=GW_TOLQ0) then 2134 write(msg,'(5a)')& 2135 'The list of q-points to be calculated contains the Gamma point, ',ch10,& 2136 'however Gamma is not the first point in the list. ' ,ch10,& 2137 'Please, change your input file accordingly. ' 2138 ABI_ERROR(msg) 2139 end if 2140 2141 ! === Initialize the band structure datatype === 2142 ! * Copy KSS energies and occupations up to Ep%nbnds==Dtset%nband(:) 2143 ! TODO Recheck symmorphy and inversion 2144 bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)) 2145 2146 ABI_MALLOC(doccde,(bantot)) 2147 ABI_MALLOC(eigen,(bantot)) 2148 ABI_MALLOC(occfact,(bantot)) 2149 doccde(:)=zero; eigen(:)=zero; occfact(:)=zero 2150 2151 jj=0; ibtot=0 2152 do isppol=1,Dtset%nsppol 2153 do ikibz=1,Dtset%nkpt 2154 do ib=1,Hdr_wfk%nband(ikibz+Dtset%nkpt*(isppol-1)) 2155 ibtot=ibtot+1 2156 if (ib<=Ep%nbnds) then 2157 jj=jj+1 2158 occfact(jj)=Hdr_wfk%occ(ibtot) 2159 eigen (jj)=energies_p(ib,ikibz,isppol) 2160 end if 2161 end do 2162 end do 2163 end do 2164 ABI_FREE(energies_p) 2165 2166 ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 2167 ! the symmetry operations in the old GW code (symmorphy and inversion) 2168 ltest = (ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz)) < tol6)) 2169 if (.not. ltest) then 2170 do jj=1,Kmesh%nibz 2171 write(std_out, *)"wtk dtset vs kmesh:", dtset%wtk(jj), kmesh%wt(jj) 2172 end do 2173 end if 2174 ABI_CHECK(ltest, 'Mismatch between Dtset%wtk and Kmesh%wt') 2175 2176 ABI_MALLOC(npwarr,(Hdr_wfk%nkpt)) 2177 npwarr(:)=Ep%npwwfn 2178 2179 call ebands_init(bantot,ks_ebands,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& 2180 & doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 2181 & Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 2182 & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 2183 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 2184 2185 ! TODO modify outkss in order to calculate the eigenvalues also if NSCF calculation. 2186 ! this fails simply because in case of NSCF occ are zero 2187 !ltest=(ALL(ABS(occfact-ks_ebands%occ)<1.d-2)) 2188 !call assert(ltest,'difference in occfact') 2189 !write(std_out,*)MAXVAL(ABS(occfact(:)-ks_ebands%occ(:))) 2190 2191 !TODO call ebands_update_occ here 2192 !$call ebands_update_occ(ks_ebands,spinmagntarget,Dtset%prtvol) 2193 2194 ABI_FREE(doccde) 2195 ABI_FREE(eigen) 2196 ABI_FREE(npwarr) 2197 2198 ! Initialize abinit header for the screening part 2199 call hdr_init(ks_ebands,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl) 2200 2201 ! Get Pawrhoij from the header. 2202 ABI_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw)) 2203 if (Dtset%usepaw==1) then 2204 call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 2205 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 2206 end if 2207 call Hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 2208 2209 ABI_FREE(occfact) 2210 call pawrhoij_free(Pawrhoij) 2211 ABI_FREE(Pawrhoij) 2212 2213 ! ==== Setup of extrapolar technique ==== 2214 Ep%gwcomp = Dtset%gwcomp; Ep%gwencomp = Dtset%gwencomp 2215 if (Ep%gwcomp == 1) then 2216 write(msg,'(a,f8.2,a)')' Using the completeness correction with gwencomp ',Ep%gwencomp*Ha_eV,' [eV] ' 2217 call wrtout(std_out, msg) 2218 end if 2219 2220 ! Final compatibility tests 2221 if (ANY(ks_ebands%istwfk /= 1)) then 2222 ABI_WARNING('istwfk/=1 is still under development') 2223 end if 2224 2225 ltest = (ks_ebands%mband == Ep%nbnds .and. ALL(ks_ebands%nband == Ep%nbnds)) 2226 ABI_CHECK(ltest,'BUG in definition of ks_ebands%nband') 2227 2228 if (Ep%gwcomp==1 .and. Ep%spmeth>0) then 2229 ABI_ERROR("Hilbert transform and extrapolar method are not compatible") 2230 end if 2231 2232 DBG_EXIT('COLL') 2233 2234 end subroutine setup_screening