TABLE OF CONTENTS
ABINIT/m_gwr_driver [ Modules ]
NAME
m_gwr_driver
FUNCTION
Driver for GWR calculations
COPYRIGHT
Copyright (C) 2021-2024 ABINIT group (MG) This file is distributed under the terms of the APACHE license version 2.0, see ~abinit/COPYING or https://www.apache.org/licenses/LICENSE-2.0 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_gwr_driver 23 24 #ifdef HAVE_MPI2 25 use mpi 26 #endif 27 use defs_basis 28 use defs_wvltypes 29 use m_errors 30 use m_abicore 31 use m_xmpi 32 use m_xomp 33 use m_hdr 34 use libxc_functionals 35 use m_crystal 36 use m_ebands 37 use m_dtset 38 use m_dtfil 39 use m_wfk 40 use m_distribfft 41 use m_nctk 42 use, intrinsic :: iso_c_binding 43 44 use defs_datatypes, only : pseudopotential_type, ebands_t 45 use defs_abitypes, only : MPI_type 46 use m_time, only : timab 47 use m_io_tools, only : file_exists, open_file, get_unit, iomode_from_fname 48 use m_time, only : cwtime, cwtime_report, sec2str 49 use m_fstrings, only : strcat, sjoin, ftoa, itoa, string_in 50 use m_slk, only : matrix_scalapack 51 use m_fftcore, only : print_ngfft, get_kg 52 use m_fft, only : fourdp 53 use m_ioarr, only : read_rhor 54 use m_energies, only : energies_type, energies_init 55 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 56 use m_pawang, only : pawang_type 57 use m_pawrad, only : pawrad_type 58 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 59 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 60 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 61 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print 62 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, & 63 pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack 64 use m_pawdij, only : pawdij, symdij_all 65 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 66 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 67 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free, paw_rho_tw_g 68 use m_kg, only : getph !, getcut 69 use m_wfd, only : test_charge 70 use m_pspini, only : pspini 71 use m_paw_correlations,only : pawpuxinit 72 use m_paw_dmft, only : paw_dmft_type 73 use m_paw_sphharm, only : setsym_ylm 74 use m_paw_mkrho, only : denfgr 75 use m_paw_nhat, only : nhatgrid, pawmknhat 76 use m_paw_tools, only : chkpawovlp, pawprt 77 use m_paw_denpot, only : pawdenpot 78 use m_paw_init, only : pawinit, paw_gencond 79 use m_pawcprj, only : pawcprj_type, pawcprj_free, pawcprj_alloc ! , paw_overlap 80 use m_ksdiago, only : ugb_t 81 !use m_fock, only : fock_type, fock_init, fock_destroy, fock_ACE_destroy, fock_common_destroy, & 82 ! fock_BZ_destroy, fock_update_exc, fock_updatecwaveocc 83 use m_mkrho, only : prtrhomxmn 84 use m_melemts, only : melflags_t 85 use m_setvtr, only : setvtr 86 use m_vhxc_me, only : calc_vhxc_me 87 use m_gwr, only : gwr_t 88 !use m_ephtk, only : ephtk_update_ebands 89 90 implicit none 91 92 private
m_gwr_driver/cc4s_gamma [ Functions ]
[ Top ] [ m_gwr_driver ] [ Functions ]
NAME
cc4s_gamma
FUNCTION
Interface with CC4S code.
Compute <b1,k|e^{-iGr}|b2,k> matrix elements and store them to disk
INPUTS
m_gwr_driver/cc4s_write_eigens [ Functions ]
[ Top ] [ m_gwr_driver ] [ Functions ]
NAME
cc4s_write_eigens
FUNCTION
Write eigenvalues in CC4S format. Only master proc should call this routine.
INPUTS
m_gwr_driver/gwr_driver [ Functions ]
[ Top ] [ m_gwr_driver ] [ Functions ]
NAME
gwr_driver
FUNCTION
Main routine for GWR calculations.
INPUTS
acell(3)=Length scales of primitive translations (bohr) codvsn=Code version dtfil<datafiles_type>=Variables related to files. dtset<dataset_type>=All input variables for this dataset. pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. Before entering the first time in the routine, 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 bethe_salpeter, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. xred(3,natom)=Reduced atomic coordinates.
NOTES
ON THE 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
147 subroutine gwr_driver(codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, xred) 148 149 !Arguments ------------------------------------ 150 !scalars 151 character(len=8),intent(in) :: codvsn 152 type(datafiles_type),intent(in) :: dtfil 153 type(dataset_type),intent(inout) :: dtset 154 type(pawang_type),intent(inout) :: pawang 155 type(pseudopotential_type),intent(inout) :: psps 156 !arrays 157 real(dp),intent(in) :: xred(3,dtset%natom) 158 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 159 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 160 161 !Local variables ------------------------------ 162 !scalars 163 integer,parameter :: master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0, optrhoij1 = 1 164 integer :: ii, comm, nprocs, my_rank, mgfftf, nfftf, omp_ncpus, work_size, nks_per_proc 165 integer :: ierr, spin, ik_ibz, nband_k, iomode__, color, io_comm 166 real(dp) :: eff, mempercpu_mb, max_wfsmem_mb, nonscal_mem 167 real(dp) :: ecore, ecut_eff, ecutdg_eff, cpu, wall, gflops, diago_cpu, diago_wall, diago_gflops 168 logical, parameter :: is_dfpt = .false. 169 logical :: read_wfk, write_wfk, cc4s_task, rectangular, rdm_update, call_pawinit 170 character(len=500) :: msg 171 character(len=fnlen) :: wfk_path, den_path, kden_path, out_path 172 type(hdr_type) :: wfk_hdr, den_hdr, kden_hdr, owfk_hdr 173 type(crystal_t) :: cryst, den_cryst, wfk_cryst 174 type(ebands_t) :: ks_ebands, owfk_ebands 175 type(pawfgr_type) :: pawfgr 176 type(wvl_data) :: wvl 177 type(mpi_type) :: mpi_enreg_seq 178 type(gwr_t) :: gwr 179 type(wfk_t) :: owfk 180 !arrays 181 real(dp), parameter :: k0(3) = zero 182 integer :: cplex, cplex_dij, cplex_rhoij 183 integer :: gnt_option,has_dijU,has_dijso,ider,izero 184 integer :: istep, moved_atm_inside, moved_rhor, n3xccc, sc_mode 185 !integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene 186 integer :: ndij !,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot 187 integer :: ngrvdw, nhatgrdim, nkxc, nkxc1, nspden_rhoij, optene, nzlmopt 188 integer :: optcut, optgr0, optgr1, optgr2, optrad, psp_gencond, option 189 integer :: rhoxsp_method, usexcnhat !, use_umklp 190 real(dp) :: compch_fft, compch_sph !,r_s,rhoav,alpha 191 !real(dp) :: drude_plsmf !,my_plsmf,ecut_eff,ecutdg_eff,ehartree 192 real(dp) :: gsqcutc_eff, gsqcutf_eff, gsqcut_shp 193 real(dp) :: vxcavg !,vxcavg_qp ucvol, 194 real(dp) :: gw_gsq !, gsqcut, gwc_gsq, gwx_gsq, 195 type(energies_type) :: KS_energies 196 type(melflags_t) :: KS_mflags 197 type(paw_dmft_type) :: Paw_dmft 198 type(ugb_t) :: ugb 199 type(xmpi_pool2d_t) :: diago_pool 200 !type(fock_type),pointer :: fock => null() 201 !arrays 202 integer :: ngfftc(18),ngfftf(18),units(2) 203 integer,allocatable :: nq_spl(:), l_size_atm(:) 204 integer,allocatable :: tmp_kstab(:,:,:), npwarr_ik(:), gvec_(:,:), istwfk_ik(:), nband_iks(:,:) 205 real(dp) :: strsxc(6), diago_info(3, dtset%nkpt, dtset%nsppol),tsec(2) 206 real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:) 207 real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:) 208 real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:), ks_vtrial(:,:), ks_vxc(:,:) 209 real(dp),allocatable :: ks_taur(:,:) !, ks_vxctau(:,:), xcctau3d(:) 210 real(dp),allocatable :: kxc(:,:), ph1d(:,:), ph1df(:,:) !qp_kxc(:,:), 211 real(dp),allocatable :: vpsp(:), xccc3d(:), dijexc_core(:,:,:) !, dij_hf(:,:,:) 212 real(dp),allocatable :: eig_k(:), occ_k(:) 213 real(dp),contiguous,pointer :: cg_k_ptr(:,:) 214 type(Paw_an_type),allocatable :: KS_paw_an(:) 215 type(Paw_ij_type),allocatable :: KS_paw_ij(:) 216 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 217 type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:) 218 type(pawpwff_t),allocatable :: Paw_pwff(:) 219 !type(pawcprj_type),allocatable :: cprj_k(:,:) 220 221 !************************************************************************ 222 223 ! This part performs the initialization of the basic objects used to perform e-ph calculations: 224 ! 225 ! 1) Crystal structure `cryst` 226 ! 2) Ground state band energies: `ks_ebands` 227 ! 5) Pseudos and PAW basic objects. 228 ! 229 ! Once we have these objects, we can call specialized routines for e-ph calculations. 230 ! Notes: 231 ! 232 ! * Any modification to the basic objects mentioned above should be done here (e.g. change of efermi) 233 ! * This routines shall not allocate big chunks of memory. The CPU-demanding sections should be 234 ! performed in the subdriver that will employ different MPI distribution schemes optimized for that particular task. 235 236 ! abirules! 237 if (.False.) write(std_out,*)xred 238 239 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 240 units(:) = [std_out, ab_out] 241 242 call cwtime(cpu, wall, gflops, "start") 243 244 ! write(msg,'(7a)')& 245 ! ' SIGMA: Calculation of the GW corrections ',ch10,ch10,& 246 ! ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,& 247 ! ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.' 248 ! call wrtout(units, msg) 249 ! 250 #if defined HAVE_GW_DPC 251 write(msg,'(a,i2,a)')'.Using double precision arithmetic; gwpc = ',gwpc,ch10 252 #else 253 write(msg,'(a,i2,a)')'.Using single precision arithmetic; gwpc = ',gwpc,ch10 254 #endif 255 call wrtout(units, msg) 256 257 ! autoparal section 258 ! TODO: This just to activate autoparal in AbiPy. Lot of things should be improved. 259 if (dtset%max_ncpus /= 0) then 260 write(ab_out,'(a)')"--- !Autoparal" 261 write(ab_out,"(a)")"# Autoparal section for GWR runs" 262 write(ab_out,"(a)") "info:" 263 write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal 264 write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 265 write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt 266 write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol 267 write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor 268 write(ab_out,"(a,i0)")" mband: ",dtset%mband 269 write(ab_out,"(3a)") " gwr_task: '",trim(dtset%gwr_task),"'" 270 271 if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then 272 work_size = dtset%nkpt * dtset%nsppol * dtset%mpw 273 max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp 274 else 275 work_size = dtset%gwr_ntau * dtset%nkpt * dtset%nsppol * dtset%mpw 276 max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp 277 end if 278 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI. 279 nonscal_mem = zero 280 281 ! List of configurations. 282 ! Assuming an OpenMP implementation with perfect speedup! 283 write(ab_out,"(a)")"configurations:" 284 285 do ii=1,dtset%max_ncpus 286 nks_per_proc = work_size / ii 287 nks_per_proc = nks_per_proc + mod(work_size, ii) 288 eff = (one * work_size) / (ii * nks_per_proc) 289 ! Add the non-scalable part and increase by 10% to account for other datastructures. 290 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 291 do omp_ncpus=1,1 !xomp_get_max_threads() 292 write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 293 write(ab_out,"(a,i0)")" mpi_ncpus: ",ii 294 write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus 295 write(ab_out,"(a,f12.9)")" efficiency: ",eff 296 write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 297 end do 298 end do 299 write(ab_out,'(a)')"..." 300 ABI_STOP("Stopping now!") 301 end if 302 303 cryst = dtset%get_crystal(img=1) 304 305 ! Some variables need to be initialized/nullify at start 306 usexcnhat = 0 307 call energies_init(KS_energies) 308 309 den_path = dtfil%fildensin; wfk_path = dtfil%fnamewffk; kden_path = dtfil%filkdensin 310 311 if (my_rank == master) then 312 ! Initialize filenames. Accept files in Fortran or in netcdf format. 313 if (nctk_try_fort_or_ncfile(den_path, msg) /= 0) then 314 ABI_ERROR(sjoin("Cannot find DEN file:", den_path, ". Error:", msg)) 315 end if 316 call wrtout(units, sjoin("- Reading GS density from: ", den_path)) 317 318 if (dtset%usekden == 1) then 319 if (nctk_try_fort_or_ncfile(kden_path, msg) /= 0) then 320 ABI_ERROR(sjoin("Cannot find KDEN file:", kden_path, ". Error:", msg)) 321 end if 322 call wrtout(units, sjoin("- Reading KDEN kinetic energy density from: ", kden_path)) 323 end if 324 call wrtout(ab_out, ch10//ch10) 325 end if ! master 326 327 ! Broadcast filenames (needed if we are using netcdf files) 328 call xmpi_bcast(den_path, master, comm, ierr) 329 call xmpi_bcast(kden_path, master, comm, ierr) 330 331 ! TODO: FFT meshes for DEN/POT should be initialized from the DEN file instead of the dtset. 332 ! Interpolating the DEN indeed breaks degeneracies in the vxc matrix elements. 333 334 call pawfgr_init(pawfgr, dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, & 335 gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=cryst%gmet, k0=k0) 336 337 call print_ngfft(ngfftc, header='Coarse FFT mesh for the wavefunctions') 338 call print_ngfft(ngfftf, header='Dense FFT mesh for densities and potentials') 339 340 ! Fake MPI_type for the sequential part. 341 call initmpi_seq(mpi_enreg_seq) 342 call init_distribfft_seq(mpi_enreg_seq%distribfft, 'c', ngfftc(2), ngfftc(3), 'all') 343 call init_distribfft_seq(mpi_enreg_seq%distribfft, 'f', ngfftf(2), ngfftf(3), 'all') 344 345 ! =========================================== 346 ! === Open and read pseudopotential files === 347 ! =========================================== 348 call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm) 349 350 ! ============================ 351 ! ==== PAW initialization ==== 352 ! ============================ 353 if (dtset%usepaw == 1) then 354 call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred) 355 356 cplex_dij = dtset%nspinor; cplex = 1; ndij = 1 357 358 ABI_MALLOC(ks_pawrhoij, (cryst%natom)) 359 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, & 360 nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=dtset%pawcpxocc) 361 call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab) 362 363 ! Test if we have to call pawinit 364 gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2 365 call paw_gencond(dtset, gnt_option, "test", call_pawinit) 366 !call_pawinit = .True. 367 368 if (psp_gencond == 1 .or. call_pawinit) then 369 call timab(553, 1, tsec) 370 gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2 371 call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, & 372 psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, & 373 dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero) 374 call timab(553,2,tsec) 375 376 ! Update internal values 377 call paw_gencond(dtset, gnt_option, "save", call_pawinit) 378 else 379 if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij = 2 380 if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2 381 end if 382 383 psps%n1xccc = maxval(pawtab(1:cryst%ntypat)%usetcore) 384 385 ! Initialize optional flags in Pawtab to zero 386 ! Cannot be done in Pawinit since the routine is called only if some parts. are changed 387 pawtab(:)%has_nabla = 0 388 pawtab(:)%lamb_shielding = zero 389 390 call setsym_ylm(cryst%gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot) 391 392 ! Initialize and compute data for DFT+U 393 Paw_dmft%use_dmft = Dtset%usedmft 394 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla, & 395 is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol, & 396 Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa) 397 398 if (my_rank == master) call pawtab_print(Pawtab) 399 400 ! Get Pawrhoij from the header of the WFK file. 401 !call pawrhoij_copy(wfk_hdr%pawrhoij, KS_Pawrhoij) 402 403 ! Evaluate form factor of radial part of phi.phj-tphi.tphj. 404 gw_gsq = max(Dtset%ecutsigx, Dtset%ecuteps) / (two*pi**2) 405 406 ! Set up q-grid, make qmax 20% larger than largest expected. 407 ABI_MALLOC(nq_spl, (Psps%ntypat)) 408 ABI_MALLOC(qmax, (Psps%ntypat)) 409 qmax = SQRT(gw_gsq)*1.2d0 410 nq_spl = Psps%mqgrid_ff 411 ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax 412 413 rhoxsp_method = 1 ! Arnaud-Alouani (default in sigma) 414 !rhoxsp_method = 2 ! Shiskin-Kresse 415 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 416 417 ABI_MALLOC(paw_pwff, (psps%ntypat)) 418 call pawpwff_init(Paw_pwff, rhoxsp_method, nq_spl, qmax, cryst%gmet, pawrad, pawtab, psps) 419 420 ABI_FREE(nq_spl) 421 ABI_FREE(qmax) 422 423 ! Variables/arrays related to the fine FFT grid 424 ABI_CALLOC(ks_nhat, (nfftf, Dtset%nspden)) 425 426 ABI_MALLOC(pawfgrtab, (cryst%natom)) 427 call pawtab_get_lsize(pawtab, l_size_atm, cryst%natom, cryst%typat) 428 429 cplex = 1 430 call pawfgrtab_init(pawfgrtab, cplex, l_size_atm, dtset%nspden, dtset%typat) 431 ABI_FREE(l_size_atm) 432 compch_fft=greatest_real 433 usexcnhat = maxval(Pawtab(:)%usexcnhat) 434 ! * 0 if Vloc in atomic data is Vbare (Blochl's formulation) 435 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation) 436 call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat))) 437 ! 438 ! Identify parts of the rectangular grid where the density has to be calculated 439 optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm 440 if (Dtset%pawcross==1) optrad=1 441 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 442 443 call nhatgrid(cryst%atindx1, cryst%gmet, cryst%natom, cryst%natom, cryst%nattyp, ngfftf, cryst%ntypat,& 444 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 445 446 call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol) 447 448 else 449 ABI_MALLOC(Paw_pwff, (0)) 450 ABI_MALLOC(Pawfgrtab, (0)) 451 end if ! End of PAW Initialization 452 453 ! Allocate these arrays anyway, since they are passed to subroutines. 454 ABI_MALLOC_IFNOT(ks_nhat, (nfftf, 0)) 455 ABI_MALLOC_IFNOT(dijexc_core, (1, 1, 0)) 456 457 !============================================= 458 ! Read density and compare crystal structures 459 ! ============================================ 460 ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden)) 461 462 call read_rhor(den_path, cplex1, dtset%nspden, nfftf, ngfftf, dtset%usepaw, mpi_enreg_seq, ks_rhor, & 463 den_hdr, ks_pawrhoij, comm, allow_interp=.False.) 464 465 den_cryst = den_hdr%get_crystal() 466 if (cryst%compare(den_cryst, header=" Comparing input crystal with DEN crystal") /= 0) then 467 ABI_ERROR("Crystal structure from input and from DEN file do not agree! Check messages above!") 468 end if 469 call den_cryst%free(); call den_hdr%free() 470 471 ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden)) 472 if (dtset%usekden == 1) then 473 call read_rhor(kden_path, cplex1, dtset%nspden, nfftf, ngfftf, 0, mpi_enreg_seq, ks_taur, & 474 kden_hdr, ks_pawrhoij, comm, allow_interp=.False.) 475 call kden_hdr%free() 476 call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_taur, optrhor=1, ucvol=cryst%ucvol) 477 end if 478 479 !======================================== 480 !==== Additional computation for PAW ==== 481 !======================================== 482 nhatgrdim = 0 483 if (dtset%usepaw == 1) then 484 ! Calculate the compensation charge nhat. 485 if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc 486 cplex = 1; ider = 2 * nhatgrdim; izero = 0 487 if (nhatgrdim > 0) then 488 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim)) 489 end if 490 if (nhatgrdim == 0) then 491 ABI_MALLOC(ks_nhatgr,(0,0,0)) 492 end if 493 494 call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,& 495 Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 496 Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,& 497 Cryst%ucvol,dtset%usewvl,Cryst%xred) 498 499 ! === Evaluate onsite energies, potentials, densities === 500 ! * Initialize variables/arrays related to the PAW spheres. 501 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 502 ABI_MALLOC(KS_paw_ij, (Cryst%natom)) 503 has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0) 504 505 call paw_ij_nullify(KS_paw_ij) 506 call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,& 507 Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 508 has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,& 509 has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1, & 510 has_dijfock=dtset%usefock) 511 512 nkxc1 = 0 513 ABI_MALLOC(KS_paw_an, (Cryst%natom)) 514 call paw_an_nullify(KS_paw_an) 515 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 516 cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1,has_vxctau=dtset%usekden) 517 518 ! Calculate onsite vxc with and without core charge. 519 nzlmopt=-1; option=0; compch_sph=greatest_real 520 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,& 521 Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 522 Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 523 Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 524 Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,& 525 Cryst%ucvol,Psps%znuclpsp) 526 527 else 528 ABI_MALLOC(ks_nhatgr, (0, 0, 0)) 529 ABI_MALLOC(ks_paw_ij, (0)) 530 ABI_MALLOC(ks_paw_an, (0)) 531 end if ! PAW 532 533 !call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,& 534 ! Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 535 536 ! For PAW, add the compensation charge on the FFT mesh, then get rho(G). 537 ! NB: ks_nhat is already included in the density stored on file so we don't need to add it.to ks_rhor 538 !if (dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat 539 540 ! TODO: Overloaded interface with units or just change the API to accept units 541 call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol) 542 call prtrhomxmn(ab_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol) 543 544 if (dtset%usekden==1) then 545 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol) 546 call prtrhomxmn(ab_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol) 547 end if 548 549 ! FFT n(r) --> n(g) 550 ABI_MALLOC(ks_rhog, (2, nfftf)) 551 call fourdp(cplex1, ks_rhog, ks_rhor(:, 1), -1, mpi_enreg_seq, nfftf, 1, ngfftf, 0) 552 553 ! Compute structure factor phases and large sphere cutoff 554 ABI_MALLOC(ph1d, (2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom)) 555 ABI_MALLOC(ph1df, (2, 3 * (2 * mgfftf + 1) * Cryst%natom)) 556 557 call getph(cryst%atindx, cryst%natom, ngfftc(1), ngfftc(2), ngfftc(3), ph1d, cryst%xred) 558 559 if (psps%usepaw == 1 .and. pawfgr%usefinegrid == 1) then 560 call getph(cryst%atindx, cryst%natom, ngfftf(1), ngfftf(2), ngfftf(3), ph1df, cryst%xred) 561 else 562 ph1df(:,:)=ph1d(:,:) 563 end if 564 565 ! The following steps have been gathered in the setvtr routine: 566 ! - get Ewald energy and Ewald forces 567 ! - compute local ionic pseudopotential vpsp 568 ! - eventually compute 3D core electron density xccc3d 569 ! - eventually compute vxc and vhartr 570 ! - set up ks_vtrial 571 ! 572 !******************************************************************* 573 !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 574 !******************************************************************* 575 576 ngrvdw = 0 577 ABI_MALLOC(grvdw, (3, ngrvdw)) 578 ABI_MALLOC(grchempottn, (3, cryst%natom)) 579 ABI_MALLOC(grewtn, (3, cryst%natom)) 580 nkxc = 0 581 if (dtset%nspden == 1) nkxc = 2 582 if (dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!! 583 ! In case of MGGA, fxc and kxc are not available and we dont need them (for now ...) 584 if (dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0 585 if (nkxc /= 0) then 586 ABI_MALLOC(kxc, (nfftf, nkxc)) 587 end if 588 589 n3xccc = 0; if (psps%n1xccc /= 0) n3xccc = nfftf 590 ABI_MALLOC(xccc3d, (n3xccc)) 591 ABI_MALLOC(ks_vhartr, (nfftf)) 592 ABI_MALLOC(ks_vtrial, (nfftf, dtset%nspden)) 593 ABI_MALLOC(vpsp, (nfftf)) 594 ABI_MALLOC(ks_vxc, (nfftf, dtset%nspden)) 595 596 ! TODO: I don't think direct diago can be used with mega-GGA due to the functional derivative wrt KS states. 597 ! TB-BK should be OK though. 598 599 !ABI_MALLOC(ks_vxctau, (nfftf, dtset%nspden * dtset%usekden)) 600 !ABI_MALLOC(xcctau3d, (n3xccc * dtset%usekden)) 601 !ABI_FREE(ks_vxctau) 602 !ABI_FREE(xcctau3d) 603 604 optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1 605 606 call setvtr(Cryst%atindx1,Dtset,KS_energies,cryst%gmet,cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 607 istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 608 Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 609 optene,Pawang,Pawrad,KS_pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,cryst%rmet,cryst%rprimd,strsxc,& 610 Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred, & 611 taur=ks_taur) !xcctau3d=xcctau3d, vxctau=ks_vxctau) 612 613 ABI_FREE(grvdw) 614 ABI_FREE(grchempottn) 615 ABI_FREE(grewtn) 616 617 !============================ 618 !==== Compute KS PAW Dij ==== 619 !============================ 620 if (dtset%usepaw == 1) then 621 call timab(561,1,tsec) 622 623 ! Calculate the unsymmetrized Dij. 624 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,& 625 Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 626 Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,& 627 Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 628 k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,& 629 nucdipmom=Dtset%nucdipmom) 630 631 ! Symmetrize KS Dij 632 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,& 633 Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,& 634 Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 635 636 ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 637 call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab) 638 call timab(561,2,tsec) 639 end if 640 641 call cwtime_report(" prepare gwr_driver_init", cpu, wall, gflops) 642 643 if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then 644 ! ========================================== 645 ! Direct diagonalization of the Hamiltonian 646 ! ========================================== 647 ABI_MALLOC(nband_iks, (dtset%nkpt, dtset%nsppol)) 648 ABI_MALLOC(npwarr_ik, (dtset%nkpt)) 649 ABI_MALLOC(istwfk_ik, (dtset%nkpt)) 650 istwfk_ik = 1 651 652 ! Compute npw_k from ecut so that we can update the header and redefine %mpw 653 do ik_ibz=1,dtset%nkpt 654 !if (dtset%istwfk(ik_ibz) == 2) istwfk_ik(ik_ibz) = 2 ! TODO: istwkf 2 is not yet supported. 655 call get_kg(dtset%kptns(:,ik_ibz), istwfk_ik(ik_ibz), dtset%ecut, cryst%gmet, npwarr_ik(ik_ibz), gvec_) 656 ABI_FREE(gvec_) 657 end do 658 dtset%mpw = maxval(npwarr_ik) 659 660 ! CC4S does not need to output the WFK file. 661 write_wfk = string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL") 662 663 ! Use input nband or min of npwarr_ik to set the number of bands. 664 if (string_in(dtset%gwr_task, "HDIAGO, CC4S")) nband_iks(:,:) = maxval(dtset%nband) 665 if (string_in(dtset%gwr_task, "HDIAGO_FULL, CC4S_FULL")) nband_iks(:,:) = minval(npwarr_ik) 666 cc4s_task = string_in(dtset%gwr_task, "CC4S, CC4S_FULL") 667 if (cc4s_task) then 668 ABI_CHECK_IEQ(dtset%nkpt, 1, "CC4S interface does not support more than one k-point.") 669 !ABI_CHECK(dtset%nkpt == 1 .and. all(abs(dtset%kptns(:,1)) < tol12), "CC4S requires Gamma-only sampling") 670 end if 671 672 ! Build header with new npwarr and nband. 673 owfk_ebands = ebands_from_dtset(dtset, npwarr_ik, nband=nband_iks) 674 owfk_ebands%eig = zero 675 call hdr_init(owfk_ebands, codvsn, dtset, owfk_hdr, pawtab, 0, psps, wvl%descr) 676 677 ! Change the value of istwfk taken from dtset. 678 ABI_REMALLOC(owfk_hdr%istwfk, (dtset%nkpt)) 679 owfk_hdr%istwfk(:) = istwfk_ik 680 681 ! Build pools to distribute (kpt, spin). Try to get rectangular grids in each pool to improve efficiency in slk diago. 682 rectangular = .True.; if (dtset%nkpt == 1) rectangular = .False. 683 call diago_pool%from_dims(dtset%nkpt, dtset%nsppol, comm, rectangular=rectangular) 684 diago_info = zero 685 686 !if (dtset%usefock == 1) then 687 ! ! Initialize data_type fock for the calculation. See also m_scfcv_core. 688 ! ABI_CHECK(dtset%usepaw == 0, "FOCK with PAW not coded") 689 ! !call fock_from_wfk(fock, dtset, cryst, pawang, pawfgr, pawtab) 690 !end if 691 692 if (write_wfk) then 693 ! Master writes header and Fortran record markers 694 out_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) out_path = nctk_ncify(out_path) 695 iomode__ = iomode_from_fname(out_path) 696 call wrtout(std_out, sjoin(" Writing wavefunctions to file:", out_path)) 697 if (my_rank == master) then 698 call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), xmpi_comm_self, write_hdr=.True., write_frm=.True.) 699 call owfk%close() 700 end if 701 call xmpi_barrier(comm) 702 end if 703 704 do spin=1,dtset%nsppol 705 do ik_ibz=1,dtset%nkpt 706 if (.not. diago_pool%treats(ik_ibz, spin)) cycle 707 call cwtime(diago_cpu, diago_wall, diago_gflops, "start") 708 709 nband_k = nband_iks(ik_ibz, spin) 710 call ugb%from_diago(spin, istwfk_ik(ik_ibz), dtset%kptns(:,ik_ibz), dtset%ecut, nband_k, ngfftc, nfftf, & 711 dtset, pawtab, pawfgr, ks_paw_ij, cryst, psps, ks_vtrial, eig_k, diago_pool%comm%value) 712 713 call cwtime(diago_cpu, diago_wall, diago_gflops, "stop") 714 if (diago_pool%comm%me == 0) diago_info(1, ik_ibz, spin) = diago_wall 715 call cwtime(diago_cpu, diago_wall, diago_gflops, "start") 716 717 owfk_ebands%eig(1:nband_k, ik_ibz, spin) = eig_k(1:nband_k) 718 719 if (write_wfk) then 720 ! occupancies are set to zero. Client code is responsible for recomputing occ and fermie when reading this WFK. 721 ABI_CALLOC(occ_k, (nband_k)) 722 color = merge(1, 0, ugb%my_nband > 0) 723 call xmpi_comm_split(diago_pool%comm%value, color, diago_pool%comm%me, io_comm, ierr) 724 725 if (ugb%my_nband > 0) then 726 ABI_CHECK(all(shape(ugb%cg_k) == [2, ugb%npwsp, ugb%my_nband]), "Wrong shape") 727 ABI_CHECK_IEQ(ugb%npw_k, owfk_hdr%npwarr(ik_ibz), "Wronk npw_k") 728 call c_f_pointer(c_loc(ugb%cg_k), cg_k_ptr, shape=[2, ugb%npwsp * ugb%my_nband]) 729 730 ! Reopen file inside io_comm. 731 call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), io_comm, & 732 write_hdr=.False., write_frm=.False.) 733 734 !sc_mode = merge(xmpio_single, xmpio_collective, ugb%has_idle_procs) 735 sc_mode = xmpio_collective 736 call owfk%write_band_block([ugb%my_bstart, ugb%my_bstop], ik_ibz, spin, sc_mode, & 737 kg_k=ugb%kg_k, cg_k=cg_k_ptr, & 738 eig_k=owfk_ebands%eig(:, ik_ibz, spin), occ_k=occ_k) 739 call owfk%close() 740 end if 741 call xmpi_comm_free(io_comm) 742 ABI_FREE(occ_k) 743 end if 744 745 call cwtime(diago_cpu, diago_wall, diago_gflops, "stop") 746 if (diago_pool%comm%me == 0) diago_info(2:3, ik_ibz, spin) = [diago_wall, dble(diago_pool%comm%nproc)] 747 748 if (cc4s_task) call cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, owfk_ebands, psps, pawtab, paw_pwff, ugb) 749 750 ABI_FREE(eig_k) 751 call ugb%free() 752 end do ! ik_ibz 753 end do ! spin 754 755 call wrtout(std_out, " Direct diago completed by this MPI pool. Other pools might take more time if k != 0") 756 757 call xmpi_sum_master(diago_info, master, comm, ierr) 758 if (my_rank == master) then 759 do spin=1,dtset%nsppol 760 do ik_ibz=1,dtset%nkpt 761 associate (info => diago_info(:, ik_ibz, spin)) 762 write(std_out, "(2(a,i0),5a,i0)") " ik_ibz: ", ik_ibz, ", spin: ", spin, & 763 ", diago_wall: ", trim(sec2str(info(1))), ", io_wall: ", trim(sec2str(info(2))), ", nprocs: ", int(info(3)) 764 end associate 765 end do 766 end do 767 end if 768 769 ! Collect eigenvalues for the different k-points/spins 770 do spin=1,dtset%nsppol 771 do ik_ibz=1,dtset%nkpt 772 if (diago_pool%treats(ik_ibz, spin) .and. diago_pool%comm%me /= 0) owfk_ebands%eig(:, ik_ibz, spin) = zero 773 end do 774 end do 775 call xmpi_sum(owfk_ebands%eig, comm, ierr) 776 777 call ebands_update_occ(owfk_ebands, dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.False.) 778 779 if (my_rank == master) then 780 call ebands_print_gaps(owfk_ebands, ab_out, header="KS gaps after direct diagonalization") 781 call ebands_print_gaps(owfk_ebands, std_out, header="KS gaps after direct diagonalization") 782 if (cc4s_task) call cc4s_write_eigens(owfk_ebands, dtfil) 783 end if 784 785 ABI_FREE(npwarr_ik) 786 ABI_FREE(istwfk_ik) 787 ABI_FREE(nband_iks) 788 call owfk_hdr%free(); call ebands_free(owfk_ebands); call diago_pool%free() 789 790 ! Deallocate exact exchange data at the end of the calculation 791 !if (dtset%usefock == 1) then 792 ! if (fock%fock_common%use_ACE/=0) call fock_ACE_destroy(fock%fockACE) 793 ! !call fock_common_destroy(fock%fock_common) 794 ! call fock_BZ_destroy(fock%fock_BZ) 795 ! call fock_destroy(fock) 796 ! nullify(fock) 797 !end if 798 799 else 800 ! ==================================================== 801 ! === This is the real GWR stuff once all is ready === 802 ! ==================================================== 803 ABI_CHECK(dtset%usepaw == 0, "PAW in GWR not yet implemented.") 804 read_wfk = .True. 805 if (read_wfk) then 806 if (my_rank == master) then 807 if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then 808 ABI_ERROR(sjoin("Cannot find GS WFK file:", wfk_path, ". Error:", msg)) 809 end if 810 call wrtout(units, sjoin("- Reading GS states from WFK file:", wfk_path)) 811 end if 812 813 ! Broadcast filenames (needed because they might have been changed if we are using netcdf files) 814 call xmpi_bcast(wfk_path, master, comm, ierr) 815 816 ! Construct crystal and ks_ebands from the GS WFK file. 817 ks_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr) 818 call wfk_hdr%vs_dtset(dtset) 819 820 wfk_cryst = wfk_hdr%get_crystal() 821 if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then 822 ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!") 823 end if 824 !call wfk_cryst%print(header="crystal structure from WFK file") 825 call wfk_cryst%free() 826 827 ! Make sure that ef is inside the gap if semiconductor. 828 call ebands_update_occ(ks_ebands, dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.True.) 829 830 ! Here we change the GS bands (Fermi level, scissors operator ...) 831 ! All the modifications to ebands should be done here. 832 !call ephtk_update_ebands(dtset, ks_ebands, "Ground state energies") 833 end if 834 835 call gwr%init(dtset, dtfil, cryst, psps, pawtab, ks_ebands, mpi_enreg_seq, comm) 836 if (gwr%idle_proc) goto 100 837 838 !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW === 839 ! * This part is parallelized within wfd%comm since each node has all GW wavefunctions. 840 ! * Note that vH matrix elements are calculated using the true uncutted interaction. 841 842 rdm_update = dtset%gwr_task == "GAMMA_GW" 843 844 call KS_mflags%reset() 845 if (rdm_update) then 846 KS_mflags%has_hbare=1 847 KS_mflags%has_kinetic=1 848 end if 849 KS_mflags%has_vhartree=1 850 KS_mflags%has_vxc =1 851 KS_mflags%has_vxcval =1 852 if (Dtset%usepawu /= 0 ) KS_mflags%has_vu = 1 853 if (Dtset%useexexch /= 0) KS_mflags%has_lexexch = 1 854 if (Dtset%usepaw==1 .and. Dtset%gw_sigxcore == 1) KS_mflags%has_sxcore = 1 855 ! off-diagonal elements only for SC on wavefunctions. 856 KS_mflags%only_diago = 1 857 if (rdm_update) KS_mflags%only_diago = 0 858 859 ! Load wavefunctions for Sigma_nk in gwr%kcalc_wfd. 860 call gwr%load_kcalc_wfd(wfk_path, tmp_kstab) 861 862 ! Compute gwr%ks_me matrix elements. 863 if (.not. string_in(dtset%gwr_task, "RPA_ENERGY")) then 864 ! FIXME: This routine allocates (nband, nband) matrices and should be rewritten! 865 call calc_vhxc_me(gwr%kcalc_wfd, ks_mflags, gwr%ks_me, cryst, dtset, nfftf, ngfftf, & 866 ks_vtrial, ks_vhartr, ks_vxc, psps, pawtab, ks_paw_an, pawang, pawfgrtab, ks_paw_ij, dijexc_core, & 867 ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur) 868 if (my_rank == master) call gwr%ks_me%print(header="KS matrix elements", unit=std_out) 869 end if 870 871 ABI_FREE(tmp_kstab) 872 873 if (read_wfk) then 874 ! Read wavefunctions from WFK file. 875 call gwr%read_ugb_from_wfk(wfk_path) 876 else 877 ! Diagonalize H on the fly and 878 !call gwr%get_ugb_from_vtrial(ngfftf, ks_vtrial) 879 !gwr%wfk_hdr = ? 880 end if 881 882 ! Now call high-level routines depending on gwr_task. 883 select case (dtset%gwr_task) 884 case ("RPA_ENERGY") 885 call gwr%rpa_energy() 886 case ("GAMMA_GW") 887 call gwr%gamma_gw(nfftf, ngfftf, vpsp) 888 case ("CHI0") 889 call gwr%run_chi0() 890 case ("G0W0") 891 call gwr%run_g0w0() 892 case ("G0V") 893 call gwr%build_sigxme() 894 case ("EGEW", "EGW0", "G0EW") 895 call gwr%run_energy_scf() 896 case default 897 ABI_ERROR(sjoin("Invalid gwr_task:", dtset%gwr_task)) 898 end select 899 end if 900 901 !===================== 902 !==== Free memory ==== 903 !===================== 904 100 call xmpi_barrier(comm) 905 ABI_FREE(ks_nhat) 906 ABI_FREE(ks_nhatgr) 907 ABI_FREE(dijexc_core) 908 call pawfgr_destroy(pawfgr) 909 910 if (dtset%usepaw == 1) then 911 ! Deallocation for PAW. 912 call pawrhoij_free(ks_pawrhoij) 913 ABI_FREE(ks_pawrhoij) 914 call pawfgrtab_free(pawfgrtab) 915 !ABI_FREE(pawfgrtab) 916 call paw_ij_free(ks_paw_ij) 917 ABI_FREE(ks_paw_ij) 918 call paw_an_free(ks_paw_an) 919 !ABI_FREE(ks_paw_an) 920 call pawpwff_free(Paw_pwff) 921 end if 922 923 ABI_FREE(ph1d) 924 ABI_FREE(ph1df) 925 ABI_FREE(ks_rhor) 926 ABI_FREE(ks_rhog) 927 ABI_FREE(ks_taur) 928 ABI_FREE(kxc) 929 ABI_FREE(xccc3d) 930 ABI_FREE(ks_vhartr) 931 ABI_FREE(ks_vtrial) 932 ABI_FREE(vpsp) 933 ABI_FREE(ks_vxc) 934 ! PAW 935 ABI_SFREE(paw_pwff) 936 ABI_SFREE(pawfgrtab) 937 ABI_SFREE(ks_paw_an) 938 939 call cryst%free(); call wfk_hdr%free(); call ebands_free(ks_ebands); call destroy_mpi_enreg(mpi_enreg_seq) 940 call gwr%free() 941 942 end subroutine gwr_driver