Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwr_driver
2 : : !! NAME
3 : : !! m_gwr_driver
4 : : !!
5 : : !! FUNCTION
6 : : !! Driver for GWR calculations
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2021-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! APACHE license version 2.0, see ~abinit/COPYING
12 : : !! or https://www.apache.org/licenses/LICENSE-2.0 .
13 : : !!
14 : : !! SOURCE
15 : :
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 : : use, intrinsic :: iso_c_binding
25 : : USE_MPI
26 : : use defs_basis
27 : : use defs_wvltypes
28 : : use m_errors
29 : : use m_abicore
30 : : use m_xmpi
31 : : use m_xomp
32 : : use m_hdr
33 : : use libxc_functionals
34 : : use m_crystal
35 : : use m_ebands
36 : : use m_dtset
37 : : use m_dtfil
38 : : use m_wfk
39 : : use netcdf
40 : : use m_nctk
41 : :
42 : : use defs_datatypes, only : pseudopotential_type
43 : : use defs_abitypes, only : MPI_type
44 : : use m_time, only : timab
45 : : use m_io_tools, only : file_exists, open_file, get_unit, iomode_from_fname
46 : : use m_time, only : cwtime, cwtime_report, sec2str
47 : : use m_fstrings, only : strcat, sjoin, ftoa, itoa, string_in, ltoa
48 : : use m_fftcore, only : print_ngfft, get_kg
49 : : use m_fft, only : fourdp
50 : : use m_ioarr, only : read_rhor
51 : : use m_energies, only : energies_type
52 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
53 : : use m_pawang, only : pawang_type
54 : : use m_pawrad, only : pawrad_type
55 : : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
56 : : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
57 : : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
58 : : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
59 : : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, &
60 : : pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack
61 : : use m_pawdij, only : pawdij, symdij_all
62 : : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
63 : : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
64 : : use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free, paw_rho_tw_g
65 : : use m_kg, only : getph
66 : : use m_wfd, only : wfd_t, test_charge
67 : : use m_pspini, only : pspini
68 : : use m_paw_correlations,only : pawpuxinit
69 : : use m_paw_dmft, only : paw_dmft_type
70 : : use m_paw_sphharm, only : setsym_ylm
71 : : use m_paw_mkrho, only : denfgr
72 : : use m_paw_nhat, only : nhatgrid, pawmknhat
73 : : use m_paw_tools, only : chkpawovlp, pawprt
74 : : use m_paw_denpot, only : pawdenpot
75 : : use m_paw_init, only : pawinit, paw_gencond
76 : : use m_pawcprj, only : pawcprj_type, pawcprj_free, pawcprj_alloc ! , paw_overlap
77 : : use m_pawxc, only : pawxc_get_usekden
78 : : use m_ksdiago, only : ugb_t, hyb_t
79 : : use m_mkrho, only : prtrhomxmn
80 : : use m_melemts, only : melflags_t
81 : : use m_setvtr, only : setvtr
82 : : use m_vhxc_me, only : calc_vhxc_me
83 : : use m_gwr, only : gwr_t
84 : : use m_vcoul, only : vcgen_t
85 : : use m_ephtk, only : ephtk_update_ebands
86 : : use m_pstat, only : pstat_proc
87 : :
88 : : implicit none
89 : :
90 : : private
91 : : !!***
92 : :
93 : : public :: gwr_driver
94 : :
95 : : #ifdef HAVE_MPI1
96 : : include 'mpif.h'
97 : : #endif
98 : : !!***
99 : :
100 : : contains
101 : : !!***
102 : :
103 : : !!****f* m_gwr_driver/gwr_driver
104 : : !! NAME
105 : : !! gwr_driver
106 : : !!
107 : : !! FUNCTION
108 : : !! Main routine for GWR calculations.
109 : : !!
110 : : !! INPUTS
111 : : !! acell(3)=Length scales of primitive translations (bohr)
112 : : !! codvsn=Code version
113 : : !! dtfil<datafiles_type>=Variables related to files.
114 : : !! dtset<dataset_type>=All input variables for this dataset.
115 : : !! pawang<pawang_type)>=PAW angular mesh and related data.
116 : : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
117 : : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
118 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
119 : : !! Before entering the first time in the routine, a significant part of Psps has been initialized :
120 : : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
121 : : !! and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
122 : : !! the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the
123 : : !! one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
124 : : !! xred(3,natom)=Reduced atomic coordinates.
125 : : !!
126 : : !! NOTES
127 : : !!
128 : : !! ON THE USE OF FFT GRIDS:
129 : : !! =================
130 : : !! In case of PAW:
131 : : !! ---------------
132 : : !! Two FFT grids are used:
133 : : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
134 : : !! It is defined by nfft, ngfft, mgfft, ...
135 : : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
136 : : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
137 : : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
138 : : !! In case of norm-conserving:
139 : : !! ---------------------------
140 : : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
141 : : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
142 : : !!
143 : : !! SOURCE
144 : :
145 : 0 : subroutine gwr_driver(codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, xred)
146 : :
147 : : !Arguments ------------------------------------
148 : : !scalars
149 : : character(len=8),intent(in) :: codvsn
150 : : type(datafiles_type),intent(in) :: dtfil
151 : : type(dataset_type),intent(inout) :: dtset
152 : : type(pawang_type),intent(inout) :: pawang
153 : : type(pseudopotential_type),intent(inout) :: psps
154 : : !arrays
155 : : real(dp),intent(in) :: xred(3,dtset%natom)
156 : : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
157 : : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
158 : :
159 : : !Local variables ------------------------------
160 : : !scalars
161 : : integer,parameter :: master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0, optrhoij1 = 1
162 : : integer :: ii, comm, nprocs, my_rank, mgfftf, nfftf, omp_ncpus, work_size, nks_per_proc
163 : : integer :: ierr, spin, ik_ibz, nband_k, iomode__, color, io_comm, usevxctau_paw !, kg_varid
164 : : real(dp) :: eff, mempercpu_mb, max_wfsmem_mb, nonscal_mem, el_temp
165 : : real(dp) :: ecore, ecut_eff, ecutdg_eff, cpu, wall, gflops, diago_cpu, diago_wall, diago_gflops
166 : : logical, parameter :: is_dfpt = .false.
167 : : logical :: read_wfk, write_wfk, cc4s_task, rectangular, with_pools, rdm_update, call_pawinit, cc4s_from_wfk
168 : : character(len=500) :: msg
169 : : character(len=fnlen) :: wfk_path, den_path, kden_path, out_path
170 : 0 : type(hdr_type) :: wfk_hdr, den_hdr, kden_hdr, owfk_hdr
171 [ # # # # : 0 : type(crystal_t) :: cryst, den_cryst, wfk_cryst
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
172 : 0 : type(ebands_t) :: ks_ebands, owfk_ebands, tmp_ebands
173 : : type(pawfgr_type) :: pawfgr
174 : 0 : type(wvl_data) :: wvl
175 : 0 : type(mpi_type) :: mpi_enreg_seq
176 [ # # # # : 0 : type(gwr_t) :: gwr
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
177 [ # # ]: 0 : type(wfk_t) :: owfk
178 : : !arrays
179 : : real(dp), parameter :: k0(3) = zero
180 : : integer :: cplex, cplex_dij, cplex_rhoij
181 : : integer :: gnt_option,has_dijU,has_dijso,ider,izero
182 : : integer :: istep, moved_atm_inside, moved_rhor, n3xccc, sc_mode
183 : : !integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene
184 : : integer :: ndij !,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot
185 : : integer :: ngrvdw, nhatgrdim, nkxc, nkxc1, nspden_rhoij, optene, nzlmopt
186 : : integer :: optcut, optgr0, optgr1, optgr2, optrad, psp_gencond, option
187 : : integer :: rhoxsp_method, usexcnhat !, use_umklp
188 : : real(dp) :: compch_fft, compch_sph !,r_s,rhoav,alpha
189 : : !real(dp) :: drude_plsmf !,my_plsmf,ecut_eff,ecutdg_eff,ehartree
190 : : real(dp) :: gsqcutc_eff, gsqcutf_eff, gsqcut_shp, vxcavg, gw_gsq, gs_fermie
191 : : type(energies_type) :: KS_energies
192 : : type(melflags_t) :: KS_mflags
193 : 0 : type(paw_dmft_type) :: Paw_dmft
194 [ # # # # : 0 : type(ugb_t) :: ugb
# # # # #
# ]
195 [ # # # # : 0 : type(hyb_t) :: hyb
# # # # #
# # # ]
196 : 0 : type(xmpi_pool2d_t) :: diago_pool
197 : : !arrays
198 : : integer :: ngfftc(18),ngfftf(18),units(2) !, grid3(3)
199 : 0 : integer,allocatable :: nq_spl(:), l_size_atm(:)
200 : 0 : integer,allocatable :: tmp_kstab(:,:,:), npwarr_ik(:), gvec_(:,:), istwfk_ik(:), nband_iks(:,:)
201 : 0 : real(dp) :: strsxc(6), diago_info(3, dtset%nkpt, dtset%nsppol),tsec(2)
202 : 0 : real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
203 : 0 : real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:)
204 : 0 : real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:), ks_vtrial(:,:), ks_vxc(:,:)
205 : 0 : real(dp),allocatable :: ks_taur(:,:) !, ks_vxctau(:,:), xcctau3d(:)
206 : 0 : real(dp),allocatable :: kxc(:,:), ph1d(:,:), ph1df(:,:) !qp_kxc(:,:),
207 : 0 : real(dp),allocatable :: vpsp(:), xccc3d(:), dijexc_core(:,:,:) !, dij_hf(:,:,:)
208 : 0 : real(dp),allocatable :: eig_k(:), occ_k(:)
209 : 0 : real(dp),contiguous,pointer :: cg_k_ptr(:,:)
210 : 0 : type(paw_an_type),allocatable :: KS_paw_an(:)
211 : 0 : type(paw_ij_type),allocatable :: KS_paw_ij(:)
212 : 0 : type(pawfgrtab_type),allocatable :: Pawfgrtab(:)
213 : 0 : type(pawrhoij_type),allocatable :: KS_Pawrhoij(:)
214 : 0 : type(pawpwff_t),allocatable :: Paw_pwff(:)
215 : : !type(pawcprj_type),allocatable :: cprj_k(:,:)
216 : : !************************************************************************
217 : :
218 : : ! This part performs the initialization of the basic objects used to perform e-ph calculations:
219 : : !
220 : : ! 1) Crystal structure `cryst`
221 : : ! 2) Ground state band energies: `ks_ebands`
222 : : ! 5) Pseudos and PAW basic objects.
223 : : !
224 : : ! Once we have these objects, we can call specialized routines for e-ph calculations.
225 : : ! Notes:
226 : : !
227 : : ! * Any modification to the basic objects mentioned above should be done here (e.g. change of efermi)
228 : : ! * This routines shall not allocate big chunks of memory. The CPU-demanding sections should be
229 : : ! performed in the subdriver that will employ different MPI distribution schemes optimized for that particular task.
230 : :
231 : : ! abirules!
232 : : if (.False.) write(std_out,*)xred
233 : :
234 : 0 : comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
235 [ # # ]: 0 : units(:) = [std_out, ab_out]
236 : :
237 : 0 : call cwtime(cpu, wall, gflops, "start")
238 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
239 : :
240 : : ! write(msg,'(a)')&
241 : : ! ' GWR: Calculation of the GW corrections with GWR code ',ch10,ch10,&
242 : : ! call wrtout(units, msg)
243 : : !
244 : : #if defined HAVE_GW_DPC
245 : 0 : write(msg,'(a,i2,a)')'.Using double precision arithmetic; gwpc = ',gwp,ch10
246 : : #else
247 : : write(msg,'(a,i2,a)')'.Using single precision arithmetic; gwpc = ',gwp,ch10
248 : : #endif
249 : 0 : call wrtout(units, msg)
250 : :
251 : : ! autoparal section
252 : : ! TODO: This just to activate autoparal in AbiPy. Lot of things should be improved.
253 [ # # ]: 0 : if (dtset%max_ncpus /= 0) then
254 : 0 : write(ab_out,'(a)')"--- !Autoparal"
255 : 0 : write(ab_out,"(a)")"# Autoparal section for GWR runs"
256 : 0 : write(ab_out,"(a)") "info:"
257 : 0 : write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal
258 : 0 : write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus
259 : 0 : write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt
260 : 0 : write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol
261 : 0 : write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor
262 : 0 : write(ab_out,"(a,i0)")" mband: ",dtset%mband
263 : 0 : write(ab_out,"(3a)") " gwr_task: '",trim(dtset%gwr_task),"'"
264 : :
265 [ # # ]: 0 : if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then
266 : 0 : work_size = dtset%nkpt * dtset%nsppol * dtset%mpw
267 : 0 : max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp
268 : : else
269 : 0 : work_size = dtset%gwr_ntau * dtset%nkpt * dtset%nsppol * dtset%mpw
270 : 0 : max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp
271 : : end if
272 : : ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI.
273 : 0 : nonscal_mem = zero
274 : :
275 : : ! List of configurations.
276 : : ! Assuming an OpenMP implementation with perfect speedup!
277 : 0 : write(ab_out,"(a)")"configurations:"
278 : :
279 [ # # ]: 0 : do ii=1,dtset%max_ncpus
280 : 0 : nks_per_proc = work_size / ii
281 : 0 : nks_per_proc = nks_per_proc + mod(work_size, ii)
282 : 0 : eff = (one * work_size) / (ii * nks_per_proc)
283 : : ! Add the non-scalable part and increase by 10% to account for other datastructures.
284 : 0 : mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
285 [ # # ]: 0 : do omp_ncpus=1,1 !xomp_get_max_threads()
286 : 0 : write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
287 : 0 : write(ab_out,"(a,i0)")" mpi_ncpus: ",ii
288 : 0 : write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus
289 : 0 : write(ab_out,"(a,f12.9)")" efficiency: ",eff
290 : 0 : write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
291 : : end do
292 : : end do
293 : 0 : write(ab_out,'(a)')"..."
294 : 0 : call wrtout(std_out, "Stopping now!")
295 : 0 : call abi_abort('PERS', exit_status=0, print_config=.false.)
296 : : end if
297 : :
298 [ # # # # : 0 : cryst = dtset%get_crystal(img=1)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
299 : :
300 : : ! Some variables need to be initialized/nullify at start
301 : 0 : usexcnhat = 0
302 : 0 : call KS_energies%init()
303 : :
304 : : !Get electronic temperature from dtset
305 [ # # # # : 0 : el_temp = merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
# # ]
306 : :
307 : 0 : den_path = dtfil%fildensin; wfk_path = dtfil%fnamewffk; kden_path = dtfil%filkdensin
308 : : !use_den = f (string_in(dtset%gwr_task, "CC4S_FROM_WFK")) then
309 : :
310 [ # # ]: 0 : if (my_rank == master) then
311 : : ! Initialize filenames. Accept files in Fortran or in netcdf format.
312 [ # # ]: 0 : if (nctk_try_fort_or_ncfile(den_path, msg) /= 0) then
313 : 0 : ABI_ERROR(sjoin("Cannot find DEN file:", den_path, ". Error:", msg))
314 : : end if
315 : 0 : call wrtout(units, sjoin("- Reading GS density from: ", den_path))
316 : :
317 [ # # ]: 0 : if (dtset%usekden == 1) then
318 [ # # ]: 0 : if (nctk_try_fort_or_ncfile(kden_path, msg) /= 0) then
319 : 0 : ABI_ERROR(sjoin("Cannot find KDEN file:", kden_path, ". Error:", msg))
320 : : end if
321 : 0 : call wrtout(units, sjoin("- Reading KDEN kinetic energy density from: ", kden_path))
322 : : end if
323 : 0 : call wrtout(ab_out, ch10//ch10)
324 : : end if ! master
325 : :
326 : : ! Broadcast filenames (needed if we are using netcdf files)
327 : 0 : call xmpi_bcast(den_path, master, comm, ierr)
328 : 0 : call xmpi_bcast(kden_path, master, comm, ierr)
329 : :
330 : : ! TODO: FFT meshes for DEN/POT should be initialized from the DEN file instead of the dtset.
331 : : ! Interpolating the DEN indeed breaks degeneracies in the vxc matrix elements.
332 : : call pawfgr_init(pawfgr, dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, &
333 : 0 : gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=cryst%gmet, k0=k0)
334 : :
335 [ # # ]: 0 : call print_ngfft([std_out], ngfftc, header='Coarse FFT mesh for the wavefunctions')
336 [ # # ]: 0 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh for densities and potentials')
337 : :
338 : : ! Fake MPI_type for the sequential part.
339 : 0 : call initmpi_seq(mpi_enreg_seq)
340 : 0 : call mpi_enreg_seq%distribfft%init_seq('c', ngfftc(2), ngfftc(3), 'all')
341 : 0 : call mpi_enreg_seq%distribfft%init_seq('f', ngfftf(2), ngfftf(3), 'all')
342 : :
343 : : ! ===========================================
344 : : ! === Open and read pseudopotential files ===
345 : : ! ===========================================
346 : 0 : call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm)
347 : :
348 : : ! ============================
349 : : ! ==== PAW initialization ====
350 : : ! ============================
351 [ # # ]: 0 : if (dtset%usepaw == 1) then
352 : 0 : call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred)
353 : :
354 : 0 : cplex_dij = dtset%nspinor; cplex = 1; ndij = 1
355 : :
356 [ # # # # : 0 : ABI_MALLOC(ks_pawrhoij, (cryst%natom))
# # # # #
# # # #
# ]
357 : : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, &
358 : 0 : nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=dtset%pawcpxocc)
359 : 0 : call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab)
360 : :
361 : : ! Test if we have to call pawinit
362 [ # # # # : 0 : gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2
# # ]
363 : 0 : call paw_gencond(dtset, gnt_option, "test", call_pawinit)
364 : : !call_pawinit = .True.
365 : :
366 [ # # # # ]: 0 : if (psp_gencond == 1 .or. call_pawinit) then
367 : 0 : call timab(553, 1, tsec)
368 : 0 : gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2
369 : : call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, &
370 : : psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, &
371 : 0 : dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero)
372 : 0 : call timab(553,2,tsec)
373 : :
374 : : ! Update internal values
375 : 0 : call paw_gencond(dtset, gnt_option, "save", call_pawinit)
376 : : else
377 [ # # # # ]: 0 : if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij = 2
378 [ # # # # ]: 0 : if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2
379 : : end if
380 : :
381 [ # # ]: 0 : psps%n1xccc = maxval(pawtab(1:cryst%ntypat)%usetcore)
382 : :
383 : : ! Initialize optional flags in Pawtab to zero
384 : : ! Cannot be done in Pawinit since the routine is called only if some parts. are changed
385 [ # # ]: 0 : pawtab(:)%has_nabla = 0
386 [ # # ]: 0 : pawtab(:)%lamb_shielding = zero
387 : :
388 : 0 : call setsym_ylm(cryst%gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot)
389 : :
390 : : ! Initialize and compute data for DFT+U
391 : 0 : Paw_dmft%use_dmft = Dtset%usedmft
392 : : call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla, &
393 : : is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol, &
394 : 0 : Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
395 : :
396 [ # # ]: 0 : if (my_rank == master) call pawtab_print(Pawtab)
397 : :
398 : : ! Get Pawrhoij from the header of the WFK file.
399 : : !call pawrhoij_copy(wfk_hdr%pawrhoij, KS_Pawrhoij)
400 : :
401 : : ! Evaluate form factor of radial part of phi.phj-tphi.tphj.
402 : 0 : gw_gsq = max(Dtset%ecutsigx, Dtset%ecuteps) / (two*pi**2)
403 : :
404 : : ! Set up q-grid, make qmax 20% larger than largest expected.
405 [ # # # # : 0 : ABI_MALLOC(nq_spl, (Psps%ntypat))
# # # # #
# # # ]
406 [ # # # # : 0 : ABI_MALLOC(qmax, (Psps%ntypat))
# # # # #
# # # ]
407 [ # # ]: 0 : qmax = SQRT(gw_gsq)*1.2d0
408 [ # # ]: 0 : nq_spl = Psps%mqgrid_ff
409 : : ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax
410 : :
411 : 0 : rhoxsp_method = 1 ! Arnaud-Alouani (default in sigma)
412 : : !rhoxsp_method = 2 ! Shiskin-Kresse
413 [ # # ]: 0 : if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
414 : :
415 [ # # # # : 0 : ABI_MALLOC(paw_pwff, (psps%ntypat))
# # # # #
# # # # #
# # # # ]
416 : 0 : call pawpwff_init(Paw_pwff, rhoxsp_method, nq_spl, qmax, cryst%gmet, pawrad, pawtab, psps)
417 : :
418 [ # # ]: 0 : ABI_FREE(nq_spl)
419 [ # # ]: 0 : ABI_FREE(qmax)
420 : :
421 : : ! Variables/arrays related to the fine FFT grid
422 [ # # # # : 0 : ABI_CALLOC(ks_nhat, (nfftf, Dtset%nspden))
# # # # #
# # # # #
# # # # #
# ]
423 : :
424 [ # # # # : 0 : ABI_MALLOC(pawfgrtab, (cryst%natom))
# # # # #
# # # #
# ]
425 : 0 : call pawtab_get_lsize(pawtab, l_size_atm, cryst%natom, cryst%typat)
426 : :
427 : : cplex = 1
428 : 0 : call pawfgrtab_init(pawfgrtab, cplex, l_size_atm, dtset%nspden, dtset%typat)
429 [ # # ]: 0 : ABI_FREE(l_size_atm)
430 : 0 : compch_fft=greatest_real
431 [ # # ]: 0 : usexcnhat = maxval(Pawtab(:)%usexcnhat)
432 : : ! * 0 if Vloc in atomic data is Vbare (Blochl's formulation)
433 : : ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation)
434 : 0 : call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat)))
435 : : !
436 : : ! Identify parts of the rectangular grid where the density has to be calculated
437 : 0 : optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm
438 [ # # ]: 0 : if (Dtset%pawcross==1) optrad=1
439 [ # # # # ]: 0 : if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
440 : :
441 : : call nhatgrid(cryst%atindx1, cryst%gmet, cryst%natom, cryst%natom, cryst%nattyp, ngfftf, cryst%ntypat,&
442 : 0 : optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
443 : :
444 : 0 : call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol)
445 : :
446 : : else
447 [ # # ]: 0 : ABI_MALLOC(Paw_pwff, (0))
448 [ # # ]: 0 : ABI_MALLOC(Pawfgrtab, (0))
449 : : end if ! End of PAW Initialization
450 : :
451 : : ! Allocate these arrays anyway, since they are passed to subroutines.
452 [ # # # # : 0 : ABI_MALLOC_IFNOT(ks_nhat, (nfftf, 0))
# # # # #
# # # ]
453 [ # # # # : 0 : ABI_MALLOC_IFNOT(dijexc_core, (1, 1, 0))
# # ]
454 : :
455 : : !=============================================
456 : : ! Read density and compare crystal structures
457 : : ! ============================================
458 [ # # # # : 0 : ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden))
# # # # #
# # # # #
# # ]
459 : :
460 : : call read_rhor(den_path, cplex1, dtset%nspden, nfftf, ngfftf, dtset%usepaw, mpi_enreg_seq, ks_rhor, &
461 : 0 : den_hdr, ks_pawrhoij, comm, allow_interp=.False., want_varname="density")
462 : :
463 [ # # # # : 0 : den_cryst = den_hdr%get_crystal()
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
464 [ # # ]: 0 : if (cryst%compare(den_cryst, header=" Comparing input crystal with DEN crystal") /= 0) then
465 : 0 : ABI_ERROR("Crystal structure from input and from DEN file do not agree! Check messages above!")
466 : : end if
467 : : ! Get fermie from the GS calculation.
468 : : ! NB: It might understimate the real fermi level, especially if the den was computed on a shifted k-mesh
469 : : ! at present it's only used to implement pseudobands
470 : 0 : gs_fermie = den_hdr%fermie
471 : 0 : call den_cryst%free(); call den_hdr%free()
472 : :
473 [ # # # # : 0 : ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden))
# # # # #
# # # # #
# # ]
474 [ # # ]: 0 : if (dtset%usekden == 1) then
475 : : call read_rhor(kden_path, cplex1, dtset%nspden, nfftf, ngfftf, 0, mpi_enreg_seq, ks_taur, &
476 : 0 : kden_hdr, ks_pawrhoij, comm, allow_interp=.False., want_varname="kinedr")
477 : 0 : call kden_hdr%free()
478 : 0 : call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_taur, optrhor=1, ucvol=cryst%ucvol)
479 : : end if
480 : :
481 : : !========================================
482 : : !==== Additional computation for PAW ====
483 : : !========================================
484 : 0 : nhatgrdim = 0
485 [ # # ]: 0 : if (dtset%usepaw == 1) then
486 : : ! Calculate the compensation charge nhat.
487 [ # # ]: 0 : if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc
488 : 0 : cplex = 1; ider = 2 * nhatgrdim; izero = 0
489 [ # # ]: 0 : if (nhatgrdim > 0) then
490 [ # # # # : 0 : ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
# # # # #
# # # # #
# # # # #
# # # ]
491 : : end if
492 [ # # ]: 0 : if (nhatgrdim == 0) then
493 [ # # # # ]: 0 : ABI_MALLOC(ks_nhatgr,(0,0,0))
494 : : end if
495 : :
496 : : call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
497 : : Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
498 : : Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
499 : 0 : Cryst%ucvol,dtset%usewvl,Cryst%xred)
500 : :
501 : : ! === Evaluate onsite energies, potentials, densities ===
502 : : ! Initialize variables/arrays related to the PAW spheres.
503 : : ! Initialize also lmselect (index of non-zero LM-moments of densities).
504 [ # # # # : 0 : ABI_MALLOC(KS_paw_ij, (Cryst%natom))
# # # # #
# # # # #
# # # # ]
505 [ # # ]: 0 : has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0)
506 : :
507 : 0 : call paw_ij_nullify(KS_paw_ij)
508 : : call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
509 : : Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
510 : : has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
511 : : has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1, &
512 : 0 : has_dijfock=dtset%usefock)
513 : :
514 : 0 : nkxc1 = 0
515 [ # # # # : 0 : ABI_MALLOC(KS_paw_an, (Cryst%natom))
# # # # #
# # # # #
# # # # ]
516 : 0 : call paw_an_nullify(KS_paw_an)
517 : 0 : usevxctau_paw=pawxc_get_usekden(dtset%ixc)
518 : : call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
519 : 0 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1,has_vxctau=usevxctau_paw)
520 : :
521 : : ! Calculate onsite vxc with and without core charge.
522 : 0 : nzlmopt=-1; option=0; compch_sph=greatest_real
523 : : call pawdenpot(compch_sph,el_temp,Cryst%gprimd,ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
524 : : Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_energies%paw,KS_paw_ij,&
525 : : Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,&
526 : : Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,&
527 : 0 : Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
528 : :
529 : : else
530 [ # # # # ]: 0 : ABI_MALLOC(ks_nhatgr, (0, 0, 0))
531 [ # # # # ]: 0 : ABI_MALLOC(ks_paw_ij, (0))
532 [ # # # # ]: 0 : ABI_MALLOC(ks_paw_an, (0))
533 : : end if ! PAW
534 : :
535 : : !call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
536 : : ! Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
537 : :
538 : : ! For PAW, add the compensation charge on the FFT mesh, then get rho(G).
539 : : ! NB: ks_nhat is already included in the density stored on file so we don't need to add it.to ks_rhor
540 : : !if (dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat
541 : :
542 : : ! TODO: Overloaded interface with units or just change the API to accept units
543 : 0 : call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol)
544 : 0 : call prtrhomxmn(ab_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol)
545 : :
546 [ # # ]: 0 : if (dtset%usekden==1) then
547 : 0 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol)
548 : 0 : call prtrhomxmn(ab_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol)
549 : : end if
550 : :
551 : : ! FFT n(r) --> n(g)
552 [ # # # # : 0 : ABI_MALLOC(ks_rhog, (2, nfftf))
# # # # #
# # # ]
553 : 0 : call fourdp(cplex1, ks_rhog, ks_rhor(:, 1), -1, mpi_enreg_seq, nfftf, 1, ngfftf, 0)
554 : :
555 : : ! Compute structure factor phases and large sphere cutoff
556 [ # # # # : 0 : ABI_MALLOC(ph1d, (2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom))
# # # # #
# # # # #
# # ]
557 [ # # # # : 0 : ABI_MALLOC(ph1df, (2, 3 * (2 * mgfftf + 1) * Cryst%natom))
# # # # #
# # # # #
# # ]
558 : :
559 : 0 : call getph(cryst%atindx, cryst%natom, ngfftc(1), ngfftc(2), ngfftc(3), ph1d, cryst%xred)
560 : :
561 [ # # # # ]: 0 : if (psps%usepaw == 1 .and. pawfgr%usefinegrid == 1) then
562 : 0 : call getph(cryst%atindx, cryst%natom, ngfftf(1), ngfftf(2), ngfftf(3), ph1df, cryst%xred)
563 : : else
564 [ # # # # ]: 0 : ph1df(:,:)=ph1d(:,:)
565 : : end if
566 : :
567 : : ! The following steps have been gathered in the setvtr routine:
568 : : ! - get Ewald energy and Ewald forces
569 : : ! - compute local ionic pseudopotential vpsp
570 : : ! - eventually compute 3D core electron density xccc3d
571 : : ! - eventually compute vxc and vhartr
572 : : ! - set up ks_vtrial
573 : : !
574 : : !*******************************************************************
575 : : !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
576 : : !*******************************************************************
577 : :
578 : 0 : ngrvdw = 0
579 [ # # ]: 0 : ABI_MALLOC(grvdw, (3, ngrvdw))
580 [ # # # # : 0 : ABI_MALLOC(grchempottn, (3, cryst%natom))
# # # # #
# # # ]
581 [ # # # # : 0 : ABI_MALLOC(grewtn, (3, cryst%natom))
# # # # ]
582 : 0 : nkxc = 0
583 [ # # ]: 0 : if (dtset%nspden == 1) nkxc = 2
584 [ # # ]: 0 : if (dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!!
585 : : ! In case of MGGA, fxc and kxc are not available and we dont need them (for now ...)
586 [ # # # # ]: 0 : if (dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0
587 [ # # ]: 0 : if (nkxc /= 0) then
588 [ # # # # : 0 : ABI_MALLOC(kxc, (nfftf, nkxc))
# # # # #
# # # # #
# # ]
589 : : end if
590 : :
591 [ # # ]: 0 : n3xccc = 0; if (psps%n1xccc /= 0) n3xccc = nfftf
592 [ # # # # : 0 : ABI_MALLOC(xccc3d, (n3xccc))
# # # # #
# # # ]
593 [ # # # # : 0 : ABI_MALLOC(ks_vhartr, (nfftf))
# # # # #
# # # ]
594 [ # # # # : 0 : ABI_MALLOC(ks_vtrial, (nfftf, dtset%nspden))
# # # # #
# # # # #
# # ]
595 [ # # # # : 0 : ABI_MALLOC(vpsp, (nfftf))
# # # # ]
596 [ # # # # : 0 : ABI_MALLOC(ks_vxc, (nfftf, dtset%nspden))
# # # # #
# # # #
# ]
597 : :
598 : : ! TODO: I don't think direct diago can be used with mega-GGA due to the functional derivative wrt KS states.
599 : : ! TB-BK should be OK though.
600 : :
601 : : !usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
602 : : !ABI_MALLOC(ks_vxctau, (nfftf, dtset%nspden * usevxctau))
603 : : !ABI_MALLOC(xcctau3d, (n3xccc * dtset%usekden))
604 : : !ABI_FREE(ks_vxctau)
605 : : !ABI_FREE(xcctau3d)
606 : :
607 : 0 : optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1
608 : :
609 : : call setvtr(Cryst%atindx1,Dtset,KS_energies,cryst%gmet,cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
610 : : istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
611 : : Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
612 : : optene,Pawang,Pawrad,KS_pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,cryst%rmet,cryst%rprimd,strsxc,&
613 : : Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred, &
614 : 0 : taur=ks_taur) !xcctau3d=xcctau3d, vxctau=ks_vxctau)
615 : :
616 [ # # ]: 0 : ABI_FREE(grvdw)
617 [ # # ]: 0 : ABI_FREE(grchempottn)
618 [ # # ]: 0 : ABI_FREE(grewtn)
619 : :
620 : : !============================
621 : : !==== Compute KS PAW Dij ====
622 : : !============================
623 [ # # ]: 0 : if (dtset%usepaw == 1) then
624 : 0 : call timab(561,1,tsec)
625 : :
626 : : ! Calculate the unsymmetrized Dij.
627 : : call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
628 : : Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
629 : : Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
630 : : Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
631 : : k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,&
632 : : ks_vxc,Cryst%xred,Dtset%znucl,&
633 : 0 : nucdipmom=Dtset%nucdipmom,spinaxis=Dtset%spinaxis)
634 : :
635 : : ! Symmetrize KS Dij
636 : : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
637 : : Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,&
638 : 0 : Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
639 : :
640 : : ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
641 : 0 : call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
642 : 0 : call timab(561,2,tsec)
643 : : end if
644 : :
645 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
646 : 0 : call cwtime_report(" prepare gwr_driver_init", cpu, wall, gflops)
647 : :
648 [ # # ]: 0 : if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then
649 : : ! ==========================================
650 : : ! Direct diagonalization of the Hamiltonian
651 : : ! ==========================================
652 [ # # # # : 0 : ABI_MALLOC(nband_iks, (dtset%nkpt, dtset%nsppol))
# # # # #
# # # # #
# # ]
653 [ # # # # : 0 : ABI_MALLOC(npwarr_ik, (dtset%nkpt))
# # # # #
# # # ]
654 [ # # # # : 0 : ABI_MALLOC(istwfk_ik, (dtset%nkpt))
# # # # ]
655 [ # # ]: 0 : istwfk_ik = 1
656 : :
657 : : ! Compute npw_k from ecut so that we can update the header and redefine %mpw
658 [ # # ]: 0 : do ik_ibz=1,dtset%nkpt
659 : : !if (dtset%istwfk(ik_ibz) == 2) istwfk_ik(ik_ibz) = 2 ! TODO: istwkf 2 is not yet supported.
660 [ # # ]: 0 : call get_kg(dtset%kptns(:,ik_ibz), istwfk_ik(ik_ibz), dtset%ecut, cryst%gmet, npwarr_ik(ik_ibz), gvec_)
661 [ # # ]: 0 : ABI_FREE(gvec_)
662 : : end do
663 [ # # ]: 0 : dtset%mpw = maxval(npwarr_ik)
664 : :
665 : : ! CC4S does not need to output the WFK file.
666 : 0 : write_wfk = string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL")
667 : :
668 : : ! Use input nband or min of npwarr_ik to set the number of bands.
669 [ # # # # : 0 : if (string_in(dtset%gwr_task, "HDIAGO, CC4S")) nband_iks(:,:) = maxval(dtset%nband)
# # # # ]
670 [ # # # # : 0 : if (string_in(dtset%gwr_task, "HDIAGO_FULL, CC4S_FULL")) nband_iks(:,:) = minval(npwarr_ik)
# # # # ]
671 : 0 : cc4s_task = string_in(dtset%gwr_task, "CC4S, CC4S_FULL")
672 [ # # ]: 0 : if (cc4s_task) then
673 [ # # ]: 0 : ABI_CHECK_IEQ(dtset%nkpt, 1, "CC4S interface does not support more than one k-point.")
674 : : end if
675 : :
676 : : ! Build header with new npwarr and nband.
677 : 0 : call owfk_ebands%from_dtset(dtset, npwarr_ik, nband=nband_iks)
678 [ # # # # : 0 : owfk_ebands%eig = zero
# # ]
679 [ # # # # : 0 : owfk_ebands%istwfk = istwfk_ik
# # # # ]
680 : : !print *, "owfk_ebands%npwarr:", owfk_ebands%npwarr; stop
681 : 0 : call owfk_hdr%init(owfk_ebands, codvsn, dtset, pawtab, 0, psps, wvl%descr)
682 : :
683 : : ! Change the value of istwfk taken from dtset and set the Fermie level from gs_fermie.
684 [ # # # # : 0 : ABI_REMALLOC(owfk_hdr%istwfk, (dtset%nkpt))
# # # # #
# # # # #
# # ]
685 [ # # ]: 0 : owfk_hdr%istwfk(:) = istwfk_ik
686 : 0 : owfk_hdr%fermie = gs_fermie
687 : :
688 : : ! Build MPI pools to distribute (kpt, spin).
689 : : ! Try to get rectangular grids in each pool to improve efficiency in slk diago.
690 [ # # ]: 0 : rectangular = .True.; if (dtset%nkpt == 1) rectangular = .False.
691 : 0 : with_pools = .True.
692 [ # # # # : 0 : if (any(dtset%gwr_np_kgts /= 0) .and. (nprocs == dtset%gwr_np_kgts(2))) then
# # # # ]
693 : 0 : with_pools = .False.; call wrtout(std_out, " Deactivating pools in HDIAGO")
694 : : end if
695 : : !with_pools = .False.
696 : 0 : call diago_pool%from_dims(dtset%nkpt, dtset%nsppol, comm, with_pools, rectangular=rectangular)
697 [ # # # # : 0 : diago_info = zero
# # ]
698 : :
699 : : ! TODO: Build hyb descriptor with hybrid orbitals from WFK file.
700 : 0 : cc4s_from_wfk = .False.
701 : : if (dtset%usefock == 1 .and. cc4s_from_wfk) then
702 : : call hyb%from_wfk_file(cryst, dtfil, dtset, psps, pawtab, ngfftc, diago_pool, comm)
703 : : end if
704 : :
705 [ # # ]: 0 : if (write_wfk) then
706 : : ! Master writes the Abinit header and the Fortran record markers
707 [ # # ]: 0 : out_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) out_path = nctk_ncify(out_path)
708 : 0 : iomode__ = iomode_from_fname(out_path)
709 : 0 : call wrtout(std_out, sjoin(" Writing wavefunctions to file:", out_path))
710 [ # # ]: 0 : if (my_rank == master) then
711 : 0 : call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), xmpi_comm_self, write_hdr=.True., write_frm=.True.)
712 : 0 : call owfk%close()
713 : : end if
714 : 0 : call xmpi_barrier(comm)
715 : : end if
716 : :
717 : : ! Build H_k(g,g') for each k-point and spin and diagonalize the Hamiltonian with Scalapack/ELPA.
718 [ # # ]: 0 : do spin=1,dtset%nsppol
719 [ # # ]: 0 : do ik_ibz=1,dtset%nkpt
720 [ # # ]: 0 : if (.not. diago_pool%treats(ik_ibz, spin)) cycle
721 : 0 : nband_k = nband_iks(ik_ibz, spin)
722 : :
723 [ # # ]: 0 : if (dtset%usefock == 1 .and. .not. cc4s_from_wfk) then
724 : 0 : call wrtout(units, "Reading ugb datatype from WFK file")
725 [ # # ]: 0 : ABI_CHECK(.not. string_in(dtset%gwr_task, "CC4S_FULL"), "CC4S_FULL cannot be used with Fock, please specify nband")
726 : : call ugb%from_wfk_file(ik_ibz, spin, istwfk_ik(ik_ibz), dtset%kptns(:,ik_ibz), nband_k, dtset, &
727 [ # # ]: 0 : dtfil, cryst, eig_k, diago_pool%comm%value)
728 : :
729 : : else
730 : 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "start")
731 : : call ugb%from_diago(spin, istwfk_ik(ik_ibz), dtset%kptns(:,ik_ibz), dtset%ecut, gs_fermie, nband_k, ngfftc, nfftf, &
732 [ # # ]: 0 : dtset, pawtab, pawfgr, ks_paw_ij, cryst, psps, ks_vtrial, eig_k, hyb, diago_pool%comm%value)
733 : 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "stop")
734 : :
735 [ # # ]: 0 : if (diago_pool%comm%me == 0) diago_info(1, ik_ibz, spin) = diago_wall
736 : 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "start")
737 : : end if
738 : :
739 [ # # ]: 0 : owfk_ebands%eig(1:nband_k, ik_ibz, spin) = eig_k(1:nband_k)
740 : :
741 [ # # ]: 0 : if (write_wfk) then
742 : : ! occupancies are set to zero. Client code is responsible for recomputing occ and fermie when reading this WFK.
743 : 0 : call wrtout(std_out, sjoin(" Writing kpoint ik_ibz:", itoa(ik_ibz), ", spin: ", itoa(spin), " to disk ..."))
744 [ # # # # : 0 : ABI_CALLOC(occ_k, (nband_k))
# # # # #
# # # #
# ]
745 [ # # ]: 0 : color = merge(1, 0, ugb%my_nband > 0)
746 : 0 : call xmpi_comm_split(diago_pool%comm%value, color, diago_pool%comm%me, io_comm, ierr)
747 : : call wrtout(std_out, sjoin(" Number of MPI procs involved in IO operation:", &
748 : 0 : itoa(xmpi_comm_size(io_comm)), "/", itoa(diago_pool%comm%nproc)))
749 : :
750 [ # # ]: 0 : if (ugb%my_nband > 0) then
751 [ # # # # : 0 : ABI_CHECK(all(shape(ugb%cg_k) == [2, ugb%npwsp, ugb%my_nband]), "Wrong shape")
# # ]
752 [ # # ]: 0 : ABI_CHECK_IEQ(ugb%npw_k, owfk_hdr%npwarr(ik_ibz), "Wrong npw_k")
753 [ # # # # : 0 : call c_f_pointer(c_loc(ugb%cg_k), cg_k_ptr, shape=[2, ugb%npwsp * ugb%my_nband])
# # # # #
# # # # #
# # # # #
# ]
754 : :
755 : : ! Reopen file inside io_comm.
756 : 0 : call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), io_comm, write_hdr=.False., write_frm=.False.)
757 : : !call wrtout(std_out, "After owfk%open_write")
758 : :
759 : : ! Write my bands.
760 : : !sc_mode = merge(xmpio_single, xmpio_collective, ugb%has_idle_procs)
761 : : !sc_mode = xmpio_collective
762 : 0 : sc_mode = xmpio_single
763 : : call owfk%write_band_block([ugb%my_bstart, ugb%my_bstop], ik_ibz, spin, sc_mode, &
764 [ # # ]: 0 : kg_k=ugb%kg_k, cg_k=cg_k_ptr, eig_k=owfk_ebands%eig(:, ik_ibz, spin), occ_k=occ_k)
765 : 0 : call owfk%close()
766 : : end if
767 : 0 : call xmpi_comm_free(io_comm)
768 [ # # ]: 0 : ABI_FREE(occ_k)
769 : : end if
770 : :
771 : 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "stop")
772 [ # # # # ]: 0 : if (diago_pool%comm%me == 0) diago_info(2:3, ik_ibz, spin) = [diago_wall, dble(diago_pool%comm%nproc)]
773 : :
774 : : ! Compute and write matrix elements required by CC4S
775 [ # # ]: 0 : if (cc4s_task) call cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, owfk_ebands, psps, pawtab, paw_pwff, ugb)
776 : :
777 [ # # ]: 0 : ABI_FREE(eig_k)
778 : 0 : call ugb%free()
779 : : end do ! ik_ibz
780 : : end do ! spin
781 : 0 : call wrtout(std_out, " Direct diago completed by this MPI pool. Other pools might take more time if k != 0")
782 : :
783 : 0 : call xmpi_sum_master(diago_info, master, comm, ierr)
784 [ # # ]: 0 : if (my_rank == master) then
785 [ # # ]: 0 : do spin=1,dtset%nsppol
786 [ # # ]: 0 : do ik_ibz=1,dtset%nkpt
787 : 0 : associate (info => diago_info(:, ik_ibz, spin))
788 : 0 : write(std_out, "(2(a,i0),5a,i0)") " ik_ibz: ", ik_ibz, ", spin: ", spin, &
789 [ # # # # ]: 0 : ", diago_wall: ", trim(sec2str(info(1))), ", io_wall: ", trim(sec2str(info(2))), ", nprocs: ", int(info(3))
790 : : end associate
791 : : end do
792 : : end do
793 : : end if
794 : :
795 : : ! Collect eigenvalues for the different k-points/spins.
796 [ # # ]: 0 : do spin=1,dtset%nsppol
797 [ # # ]: 0 : do ik_ibz=1,dtset%nkpt
798 [ # # # # : 0 : if (diago_pool%treats(ik_ibz, spin) .and. diago_pool%comm%me /= 0) owfk_ebands%eig(:, ik_ibz, spin) = zero
# # ]
799 : : end do
800 : : end do
801 : 0 : call xmpi_sum(owfk_ebands%eig, comm, ierr)
802 : :
803 : 0 : call owfk_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.False.)
804 : :
805 [ # # ]: 0 : if (my_rank == master) then
806 [ # # # # ]: 0 : if (write_wfk .and. iomode__ == IO_MODE_ETSF) then
807 [ # # ]: 0 : NCF_CHECK(owfk_ebands%ncwrite_path(cryst, out_path))
808 : : end if
809 : 0 : call owfk_ebands%print_gaps(units, header="KS gaps after direct diagonalization")
810 [ # # ]: 0 : if (cc4s_task) call cc4s_write_eigens(owfk_ebands, dtfil)
811 : : end if
812 : :
813 [ # # ]: 0 : ABI_FREE(npwarr_ik)
814 [ # # ]: 0 : ABI_FREE(istwfk_ik)
815 [ # # ]: 0 : ABI_FREE(nband_iks)
816 : 0 : call owfk_hdr%free(); call owfk_ebands%free(); call hyb%free(); call diago_pool%free()
817 : :
818 [ # # ]: 0 : else if (string_in(dtset%gwr_task, "CC4S_FROM_WFK")) then
819 : : ! Read orbitals from an external WFK file and produce output files for CC4S.
820 : :
821 : : ! Construct crystal and ks_ebands from the GS WFK file.
822 [ # # # # : 0 : tmp_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
# # # # #
# # # # #
# # # # #
# # # ]
823 [ # # # # : 0 : ks_ebands = tmp_ebands%chop(1, maxval(dtset%nband))
# # # # #
# # # # #
# # # # #
# # # #
# ]
824 : :
825 : : ! Make sure that ef is inside the gap if semiconductor.
826 : : !call ks_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.True.)
827 : :
828 : 0 : call tmp_ebands%free()
829 : 0 : call wfk_hdr%vs_dtset(dtset)
830 : :
831 [ # # # # : 0 : wfk_cryst = wfk_hdr%get_crystal()
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
832 [ # # ]: 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
833 : 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
834 : : end if
835 : 0 : call wfk_cryst%free()
836 : :
837 [ # # ]: 0 : if (my_rank == master) call cc4s_write_eigens(ks_ebands, dtfil)
838 : :
839 : : ! Build MPI pools to distribute (kpt, spin).
840 : 0 : with_pools = .True.
841 [ # # # # : 0 : if (any(dtset%gwr_np_kgts /= 0) .and. (nprocs == dtset%gwr_np_kgts(2))) then
# # # # ]
842 : 0 : with_pools = .False.; call wrtout(std_out, " Deactivating pools in HDIAGO")
843 : : end if
844 : 0 : call diago_pool%from_dims(dtset%nkpt, dtset%nsppol, comm, with_pools, rectangular=.False.)
845 [ # # ]: 0 : ABI_CHECK_IEQ(dtset%nkpt, 1, "Only Gamma-point sampling is supported")
846 [ # # ]: 0 : ABI_CHECK_IEQ(dtset%nsppol, 1, "Only spin-unpolarized calculations are supported")
847 : :
848 [ # # ]: 0 : do spin=1,dtset%nsppol
849 [ # # ]: 0 : do ik_ibz=1,dtset%nkpt
850 [ # # ]: 0 : if (.not. diago_pool%treats(ik_ibz, spin)) cycle
851 : 0 : nband_k = dtset%nband(ik_ibz + (spin-1)*dtset%nkpt)
852 : : call ugb%from_wfk_file(ik_ibz, spin, dtset%istwfk(ik_ibz), dtset%kptns(:,ik_ibz), nband_k, dtset, &
853 [ # # ]: 0 : dtfil, cryst, eig_k, diago_pool%comm%value)
854 [ # # ]: 0 : ABI_FREE(eig_k)
855 : 0 : call cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, ks_ebands, psps, pawtab, paw_pwff, ugb)
856 : 0 : call ugb%free()
857 : : end do
858 : : end do
859 : 0 : call wfk_hdr%free(); call ks_ebands%free(); call diago_pool%free()
860 : :
861 : : else
862 : : ! ====================================================
863 : : ! === This is the real GWR stuff once all is ready ===
864 : : ! ====================================================
865 [ # # ]: 0 : ABI_CHECK(dtset%usepaw == 0, "PAW in GWR not yet implemented.")
866 : 0 : read_wfk = .True.
867 : : if (read_wfk) then
868 [ # # ]: 0 : if (my_rank == master) then
869 [ # # ]: 0 : if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then
870 : 0 : ABI_ERROR(sjoin("Cannot find GS WFK file:", wfk_path, ". Error:", msg))
871 : : end if
872 : 0 : call wrtout(units, sjoin("- Reading GS states from WFK file:", wfk_path))
873 : : end if
874 : :
875 : : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
876 : 0 : call xmpi_bcast(wfk_path, master, comm, ierr)
877 : :
878 : : ! Construct crystal and ks_ebands from the GS WFK file.
879 [ # # # # : 0 : tmp_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
# # # # #
# # # # #
# # # # #
# # # ]
880 [ # # # # : 0 : ks_ebands = tmp_ebands%chop(1, maxval(dtset%nband))
# # # # #
# # # # #
# # # # #
# # # #
# ]
881 : 0 : call tmp_ebands%free()
882 : 0 : call wfk_hdr%vs_dtset(dtset)
883 : :
884 [ # # # # : 0 : wfk_cryst = wfk_hdr%get_crystal()
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
885 [ # # ]: 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
886 : 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
887 : : end if
888 : : !call wfk_cryst%print(header="crystal structure from WFK file")
889 : 0 : call wfk_cryst%free()
890 : :
891 : : ! Make sure that ef is inside the gap if semiconductor.
892 : 0 : call ks_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.True.)
893 : :
894 : : ! Here we change the GS bands (Fermi level, scissors operator ...)
895 : : ! All the modifications to ebands should be done here.
896 : : !call ephtk_update_ebands(dtset, dtfil%filqpdatain, ks_ebands, "Ground state energies", comm)
897 : : end if
898 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
899 : :
900 : 0 : call gwr%init(dtset, dtfil, cryst, psps, pawtab, ks_ebands, mpi_enreg_seq, comm)
901 [ # # ]: 0 : if (gwr%idle_proc) goto 100
902 : :
903 : : !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW ===
904 : : ! * This part is parallelized within wfd%comm since each node has all GW wavefunctions.
905 : : ! * Note that vH matrix elements are calculated using the true uncutted interaction.
906 : :
907 : 0 : rdm_update = dtset%gwr_task == "GAMMA_GW"
908 : :
909 : 0 : call KS_mflags%reset()
910 [ # # ]: 0 : if (rdm_update) then
911 : 0 : KS_mflags%has_hbare=1
912 : 0 : KS_mflags%has_kinetic=1
913 : : end if
914 : 0 : KS_mflags%has_vhartree=1
915 : 0 : KS_mflags%has_vxc =1
916 : 0 : KS_mflags%has_vxcval =1
917 [ # # ]: 0 : if (Dtset%usepawu /= 0 ) KS_mflags%has_vu = 1
918 [ # # ]: 0 : if (Dtset%useexexch /= 0) KS_mflags%has_lexexch = 1
919 [ # # # # ]: 0 : if (Dtset%usepaw==1 .and. Dtset%gw_sigxcore == 1) KS_mflags%has_sxcore = 1
920 : : ! off-diagonal elements only for SC on wavefunctions.
921 : 0 : KS_mflags%only_diago = 1
922 [ # # ]: 0 : if (rdm_update) KS_mflags%only_diago = 0
923 : :
924 : : ! Load wavefunctions for Sigma_nk in gwr%kcalc_wfd.
925 [ # # ]: 0 : call gwr%load_kcalc_wfd(wfk_path, tmp_kstab)
926 : :
927 : : ! Compute gwr%ks_me matrix elements.
928 [ # # ]: 0 : if (.not. string_in(dtset%gwr_task, "RPA_ENERGY")) then
929 : : ! FIXME: This routine allocates (nband, nband) matrices and should be rewritten!
930 : : call calc_vhxc_me(gwr%kcalc_wfd, ks_mflags, gwr%ks_me, cryst, dtset, nfftf, ngfftf, &
931 : : ks_vtrial, ks_vhartr, ks_vxc, psps, pawtab, ks_paw_an, pawang, pawfgrtab, ks_paw_ij, dijexc_core, &
932 : 0 : ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur)
933 [ # # ]: 0 : if (my_rank == master) call gwr%ks_me%print(header="KS matrix elements", unit=std_out)
934 : : end if
935 : :
936 [ # # ]: 0 : ABI_FREE(tmp_kstab)
937 : :
938 : : if (read_wfk) then
939 : : ! Read wavefunctions from WFK file.
940 : 0 : call gwr%read_ugb_from_wfk(wfk_path)
941 : : else
942 : : ! Diagonalize H on the fly and
943 : : !call gwr%get_ugb_from_vtrial(ngfftf, ks_vtrial)
944 : : !gwr%wfk_hdr = ?
945 : : end if
946 : :
947 : : ! Now call high-level routines depending on gwr_task.
948 : 0 : select case (dtset%gwr_task)
949 : : case ("RPA_ENERGY")
950 : 0 : call gwr%rpa_energy()
951 : : case ("GAMMA_GW")
952 : 0 : call gwr%gamma_gw(nfftf, ngfftf, vpsp)
953 : : case ("CHI0")
954 : 0 : call gwr%run_chi0()
955 : : case ("G0W0")
956 : 0 : call gwr%run_g0w0()
957 : : case ("G0V")
958 : 0 : call gwr%build_sigxme(compute_qp=.True.)
959 : : case ("EGEW", "EGW0", "G0EW")
960 : 0 : call gwr%run_energy_scf()
961 : : case default
962 [ # # # # : 0 : ABI_ERROR(sjoin("Invalid value for gwr_task:", dtset%gwr_task))
# # # ]
963 : : end select
964 : : end if
965 : :
966 : : !=====================
967 : : !==== Free memory ====
968 : : !=====================
969 : 0 : 100 call xmpi_barrier(comm)
970 [ # # ]: 0 : ABI_FREE(ks_nhat)
971 [ # # ]: 0 : ABI_FREE(ks_nhatgr)
972 [ # # ]: 0 : ABI_FREE(dijexc_core)
973 : 0 : call pawfgr_destroy(pawfgr)
974 : :
975 [ # # ]: 0 : if (dtset%usepaw == 1) then
976 : : ! Deallocation for PAW.
977 : 0 : call pawrhoij_free(ks_pawrhoij)
978 [ # # # # : 0 : ABI_FREE(ks_pawrhoij)
# # # # #
# # # # #
# # # # ]
979 : 0 : call pawfgrtab_free(pawfgrtab)
980 : 0 : call paw_ij_free(ks_paw_ij)
981 [ # # # # : 0 : ABI_FREE(ks_paw_ij)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
982 : 0 : call paw_an_free(ks_paw_an)
983 : 0 : call pawpwff_free(Paw_pwff)
984 : : end if
985 : :
986 [ # # ]: 0 : ABI_FREE(ph1d)
987 [ # # ]: 0 : ABI_FREE(ph1df)
988 [ # # ]: 0 : ABI_FREE(ks_rhor)
989 [ # # ]: 0 : ABI_FREE(ks_rhog)
990 [ # # ]: 0 : ABI_FREE(ks_taur)
991 [ # # ]: 0 : ABI_FREE(kxc)
992 [ # # ]: 0 : ABI_FREE(xccc3d)
993 [ # # ]: 0 : ABI_FREE(ks_vhartr)
994 [ # # ]: 0 : ABI_FREE(ks_vtrial)
995 [ # # ]: 0 : ABI_FREE(vpsp)
996 [ # # ]: 0 : ABI_FREE(ks_vxc)
997 : : ! PAW stuff
998 [ # # # # : 0 : ABI_SFREE(paw_pwff)
# # # # ]
999 [ # # # # : 0 : ABI_SFREE(pawfgrtab)
# # # # #
# # # # #
# # # # #
# ]
1000 [ # # # # : 0 : ABI_SFREE(ks_paw_an)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1001 : :
1002 : 0 : call cryst%free(); call wfk_hdr%free(); call ks_ebands%free(); call destroy_mpi_enreg(mpi_enreg_seq)
1003 : : #if defined FC_NVHPC
1004 : : call wrtout(units, "- Cannot deallocate gwr datatype if FC_NVHPC, DO NOT USE DATASETS!")
1005 : : #else
1006 : 0 : call gwr%free()
1007 : : #endif
1008 : :
1009 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
1010 : :
1011 [ # # # # : 0 : end subroutine gwr_driver
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1012 : : !!***
1013 : :
1014 : : !!****f* m_gwr_driver/cc4s_write_eigens
1015 : : !! NAME
1016 : : !! cc4s_write_eigens
1017 : : !!
1018 : : !! FUNCTION
1019 : : !! Write eigenvalues in CC4S format. Only master proc should call this routine.
1020 : : !!
1021 : : !! INPUTS
1022 : :
1023 : 0 : subroutine cc4s_write_eigens(ebands, dtfil)
1024 : :
1025 : : !Arguments ------------------------------------
1026 : : type(ebands_t),intent(in) :: ebands
1027 : : type(datafiles_type),intent(in) :: dtfil
1028 : :
1029 : : !Local variables-------------------------------
1030 : : integer :: unt, ik_ibz, spin, band, gap_err
1031 : : real(dp) :: my_fermie
1032 : : character(len=500) :: msg
1033 : : character(len=fnlen) :: filepath
1034 : 0 : type(gaps_t) :: ks_gaps
1035 : : ! *************************************************************************
1036 : :
1037 : 0 : ks_gaps = ebands%get_gaps(gap_err)
1038 [ # # ]: 0 : call ks_gaps%print([std_out], header="Kohn-Sham gaps and band edges from IBZ mesh")
1039 [ # # # # : 0 : if (any(ks_gaps%ierr /= 0)) then
# # ]
1040 : 0 : ABI_ERROR("Cannot compute gaps!")
1041 : : end if
1042 : :
1043 : : ! Make sure fermi level is within the gap. I know, the case nsppol = 2 is treated in a dirty way!
1044 : 0 : my_fermie = zero
1045 [ # # ]: 0 : do spin=1,ebands%nsppol
1046 : 0 : my_fermie = ks_gaps%vb_max(spin) + (ks_gaps%cb_min(spin) - ks_gaps%vb_max(spin)) / two
1047 : : end do
1048 : 0 : my_fermie = my_fermie / ebands%nsppol
1049 : 0 : call ks_gaps%free()
1050 : :
1051 : : ! See https://manuals.cc4s.org/user-manual/objects/EigenEnergies.html
1052 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_EigenEnergies.yaml'
1053 : 0 : write(ab_out, "(3a)")ch10," Writing Eigenenergies metadata to file: ", trim(filepath)
1054 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1055 : 0 : ABI_ERROR(msg)
1056 : : end if
1057 : 0 : write(unt,'(A)') 'version: 100'
1058 : 0 : write(unt,'(A)') 'type: Tensor'
1059 : 0 : write(unt,'(A)') 'scalarType: Real64'
1060 : 0 : write(unt,'(A)') 'dimensions:'
1061 : 0 : write(unt,'(A,I0)') '- length: ',ebands%mband * ebands%nkpt * ebands%nsppol
1062 : 0 : write(unt,'(A)') ' type: State'
1063 : 0 : write(unt,'(A)') 'elements:'
1064 : 0 : write(unt,'(A)') ' type: TextFile'
1065 : 0 : write(unt,'(A)') 'unit: 1.0 # Hartree units'
1066 : 0 : write(unt,'(A)') 'metaData:'
1067 : 0 : write(unt,'(A,E22.15)') ' fermiEnergy: ',my_fermie
1068 : 0 : write(unt,'(A)') ' energies:'
1069 : :
1070 [ # # ]: 0 : do spin=1,ebands%nsppol
1071 [ # # ]: 0 : do ik_ibz=1,ebands%nkpt
1072 [ # # ]: 0 : do band=1,ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
1073 : 0 : write(unt,"(a,e22.15)") ' - ',ebands%eig(band,ik_ibz,spin)
1074 : : end do
1075 : : end do
1076 : : end do
1077 : 0 : close(unt)
1078 : :
1079 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_EigenEnergies.elements'
1080 : 0 : write(ab_out, "(3a)")ch10," Writing Eigenenergies to file: ", trim(filepath)
1081 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1082 : 0 : ABI_ERROR(msg)
1083 : : end if
1084 : :
1085 [ # # ]: 0 : do spin=1,ebands%nsppol
1086 [ # # ]: 0 : do ik_ibz=1,ebands%nkpt
1087 [ # # ]: 0 : do band=1,ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
1088 : 0 : write(unt,"(e22.15)") ebands%eig(band,ik_ibz,spin)
1089 : : end do
1090 : : end do
1091 : : end do
1092 : :
1093 : 0 : close(unt)
1094 : :
1095 [ # # # # : 0 : end subroutine cc4s_write_eigens
# # # # #
# # # # #
# # ]
1096 : : !!***
1097 : :
1098 : : !!****f* m_gwr_driver/cc4s_gamma
1099 : : !! NAME
1100 : : !! cc4s_gamma
1101 : : !!
1102 : : !! FUNCTION
1103 : : !! Interface with CC4S code.
1104 : : !! Compute <b1,k|e^{-iGr}|b2,k> matrix elements and store them to disk
1105 : : !!
1106 : : !! INPUTS
1107 : :
1108 : 0 : subroutine cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, ebands, psps, pawtab, paw_pwff, ugb)
1109 : :
1110 : : use m_numeric_tools, only : blocked_loop
1111 : : use m_gwdefs, only : GW_Q0_DEFAULT
1112 : : use m_fftcore, only : sphereboundary
1113 : : use m_fft_mesh, only : setmesh
1114 : : use m_fft, only : uplan_t
1115 : : use m_vcoul, only : vcgen_t
1116 : :
1117 : : use m_sort, only : sort_gvecs
1118 : : use m_pawpwij, only : pawpwij_t, pawpwij_init, pawpwij_free
1119 : :
1120 : : !Arguments ------------------------------------
1121 : : integer,intent(in) :: spin, ik_ibz
1122 : : type(dataset_type),intent(in) :: dtset
1123 : : type(datafiles_type),intent(in) :: dtfil
1124 : : type(crystal_t),intent(in) :: cryst
1125 : : type(ebands_t),intent(in) :: ebands
1126 : : type(Pseudopotential_type),intent(in) :: psps
1127 : : type(pawpwff_t),intent(in) :: paw_pwff(dtset%ntypat*dtset%usepaw)
1128 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1129 : : type(ugb_t),target,intent(in) :: ugb
1130 : :
1131 : : !Local variables-------------------------------
1132 : : !scalars
1133 : : integer,parameter :: mG0(3) = 0, master = 0
1134 : : integer :: nproc, my_rank, my_ib2st, npw_k, nspinor, m_npw, npwvec, ig, mpierr, fh, comm, buf_size, ierr
1135 : : integer :: band1, band1_start, batch1_size, n1dat, idat1, m_istwfk, iatom, dim_rtwg
1136 : : integer :: band2, band2_start, batch2_size, n2dat, idat2, units(2), ii, unt, nqibz_, nqbz_, nkbz_, test_unt, M_
1137 : : integer(XMPI_OFFSET_KIND) :: offset
1138 : : real(dp) :: cpu, wall, gflops, qpt(3), qbz_(3,1), gcart(3), kpt(3), max_abs_err, abs_err, my_gw_qlwl(3), mem_mb
1139 : : character(len=500) :: msg
1140 : : character(len=fnlen) :: filepath, cvx_filepath
1141 : : logical :: k_is_gamma
1142 : : logical,parameter :: trust_no_one = .False.
1143 : 0 : type(uplan_t) :: uplan_1, uplan_2, uplan_m
1144 [ # # # # : 0 : type(vcgen_t) :: vcgen
# # # # #
# ]
1145 : 0 : integer :: u_ngfft(18), u_nfft, u_mgfft, enforce_sym, method, nlmn_atm(cryst%natom)
1146 : 0 : integer,pointer :: gvec_max(:,:)
1147 : 0 : integer,allocatable,target :: m_gvec(:,:), sorted_kg_k(:,:)
1148 : 0 : complex(dp),allocatable :: ug1_batch(:,:), ur1_batch(:,:), ur2_batch(:,:), ur12_batch(:,:), ug12_batch(:,:), cwork(:)
1149 : 0 : complex(gwp),allocatable :: sqrt_vc(:), paw_rhotwg(:)
1150 : 0 : type(pawpwij_t),allocatable :: pwij(:)
1151 : 0 : type(pawcprj_type),allocatable :: cprj1(:,:)
1152 : : ! *************************************************************************
1153 : :
1154 : 0 : call cwtime(cpu, wall, gflops, "start")
1155 : :
1156 : 0 : comm = ugb%comm; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1157 [ # # ]: 0 : units = [std_out, ab_out]
1158 : 0 : npw_k = ugb%npw_k; nspinor = ugb%nspinor
1159 : :
1160 [ # # # # ]: 0 : if (dtset%prtvol > 10) call ugb%print([std_out], dtset%prtvol, header="ugb for CC4S")
1161 : :
1162 : : ! m_gvec is the g-sphere for the oscillators M computed from ecuteps (half-sphere if wavefunctions have TR).
1163 : : ! setmesh assumes g-vectors sorted by norm so use kin_sorted = True and sort ug%kg_k below.
1164 [ # # # # : 0 : kpt = dtset%kptns(:,ik_ibz); k_is_gamma = all(abs(kpt) < tol12)
# # ]
1165 [ # # ]: 0 : m_istwfk = 1; if (ugb%istwf_k == 2) m_istwfk = 2
1166 [ # # ]: 0 : call get_kg(kpt, m_istwfk, dtset%ecuteps, cryst%gmet, m_npw, m_gvec, kin_sorted=.True.)
1167 : :
1168 : : ! Setup FFT mesh
1169 [ # # ]: 0 : u_ngfft = dtset%ngfft
1170 : 0 : method = 2
1171 [ # # ]: 0 : if (dtset%fftgw==00 .or. dtset%fftgw==01) method=0
1172 [ # # ]: 0 : if (dtset%fftgw==10 .or. dtset%fftgw==11) method=1
1173 [ # # ]: 0 : if (dtset%fftgw==20 .or. dtset%fftgw==21) method=2
1174 [ # # ]: 0 : if (dtset%fftgw==30 .or. dtset%fftgw==31) method=3
1175 : 0 : enforce_sym = mod(dtset%fftgw, 10)
1176 : : ! Gamma only --> we don't need to rotate wavefunctions in the BZ
1177 [ # # ]: 0 : if (k_is_gamma) enforce_sym = 0
1178 : :
1179 [ # # ]: 0 : call sort_gvecs(npw_k, kpt, cryst%gmet, ugb%kg_k, sorted_kg_k)
1180 : :
1181 : 0 : npwvec = npw_k; gvec_max => sorted_kg_k
1182 [ # # ]: 0 : if (m_npw > npw_k) then
1183 : 0 : npwvec = m_npw; gvec_max => m_gvec
1184 : : end if
1185 [ # # ]: 0 : call setmesh(cryst%gmet, gvec_max, u_ngfft, npwvec, m_npw, npw_k, u_nfft, method, mG0, cryst, enforce_sym, unit=std_out)
1186 [ # # ]: 0 : ABI_FREE(sorted_kg_k)
1187 : :
1188 : : u_mgfft = maxval(u_ngfft(1:3))
1189 : 0 : qpt = zero
1190 : :
1191 : 0 : nqibz_ = 1; nqbz_ = 1; qbz_ = zero; nkbz_ = 1
1192 : : ! TODO: MC technique does not seem to work as expected, even in the legacy code.
1193 : : call vcgen%init(cryst, ebands%kptrlatt, nkbz_, nqibz_, nqbz_, qbz_, &
1194 : 0 : dtset%gw_rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, comm)
1195 : :
1196 : : ! NB: npweps = m_npw
1197 [ # # # # : 0 : ABI_MALLOC(sqrt_vc, (m_npw))
# # # # #
# # # ]
1198 [ # # # # ]: 0 : my_gw_qlwl(:) = GW_Q0_DEFAULT; if (dtset%gw_nqlwl > 0) my_gw_qlwl = dtset%gw_qlwl(:,1)
1199 : 0 : call vcgen%get_vc_sqrt(qpt, m_npw, m_gvec, my_gw_qlwl, cryst, sqrt_vc, comm)
1200 : : ! Override the G=0 component with the correct value set by gw_icutcoulomb method
1201 : 0 : sqrt_vc(1) = sqrt(vcgen%i_sz)
1202 : : !sqrt_vc(1) = zero
1203 : 0 : call vcgen%free()
1204 : :
1205 [ # # ]: 0 : if (my_rank == master) then
1206 : 0 : call wrtout(units, " Computing oscilator matrix elements for CC4S.")
1207 [ # # ]: 0 : call print_ngfft([std_out], u_ngfft, header='FFT mesh for wavefunctions')
1208 : :
1209 : : ! =====================
1210 : : ! Write files for CC4S
1211 : : ! =====================
1212 [ # # # # ]: 0 : if (ik_ibz == 1 .and. spin == 1) then
1213 : : ! Write g-vector files, see https://manuals.cc4s.org/user-manual/objects/GridVectors.html
1214 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_GridVectors.yaml'
1215 : 0 : write(ab_out, "(3a)")ch10," Writing Gridvectors metadata to file: ", trim(filepath)
1216 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1217 : 0 : ABI_ERROR(msg)
1218 : : end if
1219 : 0 : write(unt,'(a)') 'version: 100'
1220 : 0 : write(unt,'(a)') 'type: Tensor'
1221 : 0 : write(unt,'(a)') 'scalarType: Real64'
1222 : 0 : write(unt,'(a)') 'dimensions:'
1223 : 0 : write(unt,'(a)') ' - length: 3'
1224 : 0 : write(unt,'(a)') ' type: Vector'
1225 : 0 : write(unt,'(a,i0)') ' - length: ',m_npw
1226 : 0 : write(unt,'(a)') ' type: Momentum'
1227 : 0 : write(unt,'(a)') 'elements:'
1228 : 0 : write(unt,'(a)') ' type: TextFile'
1229 : 0 : write(unt,'(a)') 'unit: 1.0 # Bohr^-1'
1230 : : ! The last three lines correspond to the reciprocal lattice vectors (including the factor 2pi)
1231 : 0 : write(unt,'(a)') 'metaData:'
1232 : : write(unt,'(" Gi: [",E22.15,",",E22.15,",",E22.15,"]")') &
1233 : 0 : two_pi*cryst%gprimd(1,1), two_pi*cryst%gprimd(2,1), two_pi*cryst%gprimd(3,1)
1234 : : write(unt,'(" Gj: [",E22.15,",",E22.15,",",E22.15,"]")') &
1235 : 0 : two_pi*cryst%gprimd(1,2), two_pi*cryst%gprimd(2,2), two_pi*cryst%gprimd(3,2)
1236 : : write(unt,'(" Gk: [",E22.15,",",E22.15,",",E22.15,"]")') &
1237 : 0 : two_pi*cryst%gprimd(1,3), two_pi*cryst%gprimd(2,3), two_pi*cryst%gprimd(3,3)
1238 : 0 : close(unt)
1239 : :
1240 : : ! Write g-vectors (cart coords)
1241 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_GridVectors.elements'
1242 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1243 : 0 : ABI_ERROR(msg)
1244 : : end if
1245 : :
1246 : : !TODO: k-points to be implemented.
1247 [ # # ]: 0 : do ig=1,m_npw
1248 [ # # # # : 0 : gcart = two_pi * matmul(cryst%gprimd, m_gvec(:,ig))
# # # # #
# ]
1249 [ # # ]: 0 : do ii=1,3
1250 : 0 : write(unt, *) gcart(ii)
1251 : : end do
1252 : : end do
1253 : 0 : close(unt)
1254 : :
1255 : : ! https://manuals.cc4s.org/user-manual/objects/CoulombVertex.html
1256 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombVertex.yaml'
1257 : 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombVertex metadata to file: ', trim(filepath)
1258 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1259 : 0 : ABI_ERROR(msg)
1260 : : end if
1261 : 0 : write(unt,'(a)') 'version: 100'
1262 : 0 : write(unt,'(a)') 'type: Tensor'
1263 : 0 : write(unt,'(a)') 'scalarType: Complex64'
1264 : 0 : write(unt,'(a)') 'dimensions:'
1265 : 0 : write(unt,'(a,i0)') '- length: ',m_npw
1266 : 0 : write(unt,'(a)') ' type: AuxiliaryField'
1267 : 0 : write(unt,'(a,i0)') '- length: ',ugb%nband_k
1268 : 0 : write(unt,'(a)') ' type: State'
1269 : 0 : write(unt,'(a,i0)') '- length: ',ugb%nband_k
1270 : 0 : write(unt,'(a)') ' type: State'
1271 : 0 : write(unt,'(a)') 'elements:'
1272 : 0 : write(unt,'(a)') ' type: IeeeBinaryFile'
1273 : 0 : write(unt,'(a)') 'unit: 1.0 # Atomic units'
1274 : 0 : write(unt,'(a)') 'metaData:'
1275 [ # # ]: 0 : if (m_istwfk == 2) then
1276 : 0 : write(unt,'(a)') ' halfGrid: 1'
1277 : : else
1278 : 0 : write(unt,'(a)') ' halfGrid: 0'
1279 : : end if
1280 : 0 : close(unt)
1281 : :
1282 : : ! https://manuals.cc4s.org/user-manual/objects/CoulombPotential.html
1283 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombPotential.yaml'
1284 : 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombPotential metadata to file: ', trim(filepath)
1285 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1286 : 0 : ABI_ERROR(msg)
1287 : : end if
1288 : 0 : write(unt,'(a)') 'version: 100'
1289 : 0 : write(unt,'(a)') 'type: Tensor'
1290 : 0 : write(unt,'(a)') 'scalarType: Real64'
1291 : 0 : write(unt,'(a)') 'dimensions:'
1292 : 0 : write(unt,'(a,i0)') ' - length: ',m_npw
1293 : 0 : write(unt,'(a)') ' type: Momentum'
1294 : 0 : write(unt,'(a)') 'elements:'
1295 : 0 : write(unt,'(a)') ' type: TextFile'
1296 : 0 : write(unt,'(a)') 'unit: 1.0 # Atomic units '
1297 : 0 : close(unt)
1298 : :
1299 [ # # # # ]: 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombPotential.elements'
1300 : 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombPotential data to file: ', trim(filepath)
1301 [ # # ]: 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1302 : 0 : ABI_ERROR(msg)
1303 : : end if
1304 [ # # ]: 0 : do ig=1,m_npw
1305 : 0 : write(unt,*)real(sqrt_vc(ig) * conjg(sqrt_vc(ig)), kind=dp)
1306 : : !write(unt,*)real(sqrt_vc(ig)**2)
1307 : : end do
1308 : 0 : close(unt)
1309 : : end if ! ik_ibz == 1 .and. spin == 1
1310 : : end if ! my_rank == master
1311 : :
1312 : : ! Open binary file to store CoulombVertex
1313 [ # # # # ]: 0 : cvx_filepath = trim(dtfil%filnam_ds(4))//'_CoulombVertex.elements'
1314 [ # # ]: 0 : if (my_rank == master) write(ab_out, "(3a)")ch10, ' Writing CoulombVertex data to file: ', trim(cvx_filepath)
1315 : :
1316 : : #ifdef HAVE_MPI_IO
1317 : 0 : call MPI_FILE_OPEN(comm, cvx_filepath, MPI_MODE_CREATE + MPI_MODE_WRONLY, xmpio_info, fh, mpierr)
1318 : 0 : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1319 : : #else
1320 : : ABI_ERROR("CC4S interface requires MPI-IO!")
1321 : : #endif
1322 : :
1323 : : if (trust_no_one .and. my_rank == 0) then
1324 : : if (open_file("test_mg", msg, newunit=test_unt, form="formatted", status="replace", action="write") /= 0) then
1325 : : ABI_ERROR(msg)
1326 : : end if
1327 : : end if
1328 : :
1329 : : ! Define batch sizes and allocate workspace arrays.
1330 : : ! Increasing this value improves efficiency (less communication) at the price of more memory.
1331 : 0 : batch1_size = min(48, ugb%nband_k); batch2_size = min(48, ugb%nband_k) !; batch1_size = 1; batch2_size = 1
1332 : 0 : call wrtout(std_out, sjoin(" Using batch1_size:", itoa(batch1_size), ", batch2_size:", itoa(batch2_size)))
1333 : :
1334 : : mem_mb = (two * u_nfft * nspinor * batch1_size + &
1335 : : two * u_nfft * nspinor * batch2_size * two + &
1336 : 0 : two * m_npw * nspinor * batch2_size) * dp * b2Mb
1337 : 0 : call wrtout(std_out, sjoin(" Memory for workspace arrays ", ftoa(mem_mb, fmt="f8.1"), "[Mb] <<< MEM"))
1338 : :
1339 [ # # # # : 0 : ABI_MALLOC(ur1_batch, (u_nfft * nspinor, batch1_size))
# # # # #
# # # # #
# # ]
1340 [ # # # # : 0 : ABI_MALLOC(ur2_batch, (u_nfft * nspinor, batch2_size))
# # # # #
# # # # #
# # ]
1341 [ # # # # : 0 : ABI_MALLOC(ur12_batch, (u_nfft * nspinor, batch2_size))
# # # # #
# # # ]
1342 [ # # # # : 0 : ABI_MALLOC(ug12_batch, (m_npw * nspinor, batch2_size))
# # # # #
# # # # #
# # ]
1343 : :
1344 [ # # ]: 0 : if (psps%usepaw == 1) then
1345 : : ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma
1346 [ # # # # : 0 : ABI_MALLOC(pwij, (psps%ntypat))
# # # # #
# # # # #
# # # # #
# ]
1347 : 0 : call pawpwij_init(pwij, m_npw, qpt, m_gvec, cryst%rprimd, psps, pawtab, paw_pwff)
1348 [ # # ]: 0 : do iatom=1,cryst%natom
1349 : 0 : nlmn_atm(iatom) = pawtab(cryst%typat(iatom))%lmn_size
1350 : : end do
1351 [ # # # # : 0 : ABI_MALLOC(cprj1, (cryst%natom, nspinor*batch1_size))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
1352 : 0 : call pawcprj_alloc(cprj1, 0, nlmn_atm)
1353 : 0 : dim_rtwg = nspinor
1354 [ # # # # : 0 : ABI_MALLOC(paw_rhotwg, (m_npw*dim_rtwg))
# # # # #
# # # ]
1355 : : end if
1356 : :
1357 : : ! TODO:
1358 : : ! 1) take advantage of M_{b1,b2}(g) = <b1|e^{-ig.r}|b2> => M_{b1,b2}(g) = M_{b2,b1}(-g)^*
1359 : : ! once I have a better understanding of the fileformat expected by CC4S.
1360 : : ! 2) Handle parallel IO if nsppol 2 (we are inside the spin loop that is already MPI distributed!)
1361 : : ! 3) Clarify ordering of CoulombVertex (b1,b2 vs b2,b1) and eigenvalues (spin?)
1362 : : ! 4) Treatment of q--> 0 in vc_coul
1363 : : ! 4) See other TODOs below.
1364 : :
1365 : 0 : call uplan_1%init(npw_k, nspinor, batch1_size, u_ngfft, ugb%istwf_k, ugb%kg_k, dp, dtset%gpu_option)
1366 : 0 : call uplan_2%init(npw_k, nspinor, batch2_size, u_ngfft, ugb%istwf_k, ugb%kg_k, dp, dtset%gpu_option)
1367 : 0 : call uplan_m%init(m_npw, nspinor, batch2_size, u_ngfft, m_istwfk, m_gvec, dp, dtset%gpu_option)
1368 : :
1369 : 0 : M_ = m_npw
1370 : :
1371 : : ! Blocked loop over group of b1 indices. NB: Assuming bands distributed in contiguous blocks.
1372 [ # # # # : 0 : do band1_start=1, ugb%nband_k, batch1_size
# # # # ]
1373 : : ! Collect n1dat bands starting from band1_start on each proc.
1374 : 0 : n1dat = blocked_loop(band1_start, ugb%nband_k, batch1_size)
1375 : :
1376 [ # # # # ]: 0 : call ugb%mat%collect_cplx(ugb%npwsp, n1dat, [1, band1_start], ug1_batch)
1377 [ # # ]: 0 : if (psps%usepaw == 1) call ugb%collect_cprj(nspinor, n1dat, band1_start, cprj1)
1378 : :
1379 : : ! FFT: ug1_batch --> ur1_batch
1380 : 0 : call uplan_1%execute_gr(n1dat, ug1_batch(:,1), ur1_batch(:,1))
1381 [ # # # # : 0 : if (ugb%istwf_k /= 2) ur1_batch = conjg(ur1_batch) ! Not needed if k == Gamma as ur1 is real.
# # ]
1382 [ # # ]: 0 : ABI_FREE(ug1_batch)
1383 : :
1384 : : ! Blocked loop over MY group of b2 indices (contiguous blocks)
1385 [ # # # # : 0 : do band2_start=ugb%my_bstart, ugb%my_bstop, batch2_size
# # # # ]
1386 : 0 : n2dat = blocked_loop(band2_start, ugb%my_bstop, batch2_size)
1387 : 0 : my_ib2st = band2_start - ugb%my_bstart + 1
1388 : :
1389 : : ! FFT: ugb%mat --> ur2_batch for n2dat states.
1390 : 0 : call uplan_2%execute_gr(n2dat, ugb%mat%buffer_cplx(:,my_ib2st), ur2_batch(:,1))
1391 : :
1392 : : ! For each row of the submatrix, build n2dat products (band1, idat2) in r-space, then r --> g.
1393 [ # # ]: 0 : do idat1=1,n1dat
1394 : 0 : band1 = band1_start + idat1 - 1
1395 : :
1396 [ # # ]: 0 : do idat2=1,n2dat
1397 [ # # ]: 0 : ur12_batch(:,idat2) = ur1_batch(:,idat1) * ur2_batch(:,idat2)
1398 : : end do
1399 : 0 : call uplan_m%execute_rg(n2dat, ur12_batch(:,1), ug12_batch(:,1))
1400 : :
1401 [ # # ]: 0 : if (psps%usepaw == 1) then
1402 : : ! Add PAW on-site contributions
1403 [ # # ]: 0 : do idat2=1,n2dat
1404 : : associate (cprj1_kmq => cprj1(:, 1 + (idat1-1)*nspinor), &
1405 : 0 : cprj2_k => ugb%cprj_k(:, 1 + (my_ib2st+idat2-2)*nspinor)) ! NB: ugb%cprj_k(2, nspinor*my_nband)
1406 [ # # ]: 0 : paw_rhotwg = zero
1407 : 0 : call paw_rho_tw_g(cryst, pwij, m_npw, dim_rtwg, nspinor, m_gvec, cprj1_kmq, cprj2_k, paw_rhotwg)
1408 [ # # ]: 0 : ug12_batch(:,idat2) = ug12_batch(:,idat2) + paw_rhotwg
1409 : : end associate
1410 : : end do
1411 : : end if
1412 : :
1413 : : ! This to zero the matrix elements for certain (band1, band2) entries.
1414 : : ! in this case: band1 = 1, band2 = 1.
1415 : : !if (band1 == 1) then
1416 : : ! do idat2=1,n2dat
1417 : : ! ! This is the condition for band2
1418 : : ! if (band2_start + idat2 - 1 == 1) ug12_batch(:,idat2) = zero
1419 : : ! end do
1420 : : !end if
1421 : :
1422 [ # # ]: 0 : if (nspinor == 2) then
1423 : : ! Sum over spinors and repack data in the first n2dat positions to prepare IO operation.
1424 [ # # ]: 0 : do idat2=1,n2dat
1425 [ # # ]: 0 : ug12_batch(1:m_npw,idat2) = ug12_batch(1:m_npw,idat2) + ug12_batch(m_npw+1:,idat2)
1426 : : end do
1427 [ # # # # ]: 0 : do idat2=2,n2dat,2
1428 [ # # ]: 0 : ug12_batch(m_npw+1:,idat2-1) = ug12_batch(1:m_npw,idat2)
1429 : : end do
1430 : : end if
1431 : :
1432 [ # # ]: 0 : do idat2=1,n2dat
1433 : : !if (band1 == band2_start + idat2 -1) then
1434 : : ! write(std_out,*) " ug12_batch(g=0,band1,band2), band", ug12_batch(1,idat2), band1, band2_start + idat2 -1
1435 : : !end if
1436 : : ! Multiply by sqrt(vc(g))
1437 : : ! This renormalization is needed to make CC4S converge (work done by AlejandroG, FabienB, MatteoG)
1438 [ # # ]: 0 : ug12_batch(:,idat2) = ug12_batch(:,idat2) * sqrt_vc(:) / sqrt(cryst%ucvol)
1439 : : !FBru:
1440 : : ! CC4S convention:
1441 : : ! When using time-reversal, include a factor sqrt(2) for the non-zero G's
1442 : : ! Then CC4S will perform the sum only over half the grid so all the G/=0 component should be counted twice.
1443 [ # # # # ]: 0 : if (m_istwfk == 2) ug12_batch(2:,idat2) = ug12_batch(2:,idat2) * sqrt(two)
1444 : : end do
1445 : : !write(std_out,*)" max(abs(ug12_batch)):", maxval(abs(ug12_batch(:,1:n2dat)))
1446 : :
1447 : : #ifdef HAVE_MPI_IO
1448 : : ! Write ug12_batch using Stream-IO
1449 : 0 : buf_size = m_npw * n2dat
1450 : :
1451 : : ! FBru convert to longer integer before the multiplication
1452 : : offset = INT(m_npw, KIND=XMPI_OFFSET_KIND) &
1453 : : * INT( (band2_start-1) + (band1-1) * ugb%nband_k, KIND=XMPI_OFFSET_KIND) &
1454 : 0 : * INT(xmpi_bsize_dpc, KIND=XMPI_OFFSET_KIND)
1455 : :
1456 : 0 : call MPI_FILE_WRITE_AT(fh, offset, ug12_batch, buf_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1457 [ # # ]: 0 : ABI_HANDLE_MPIERR(mpierr)
1458 : :
1459 : 0 : if (my_rank == 0 .and. trust_no_one) then
1460 : : do idat2=1,n2dat
1461 : : band2 = band2_start + idat2 - 1
1462 : : write(test_unt,*)band1, band2, ug12_batch(1:M_,idat2)
1463 : : end do
1464 : : end if
1465 : : #endif
1466 : : end do ! idat1
1467 : :
1468 : : end do ! band2_start
1469 : : end do ! band1_start
1470 : :
1471 : : #ifdef HAVE_MPI_IO
1472 : 0 : call MPI_FILE_CLOSE(fh, mpierr)
1473 : 0 : ABI_CHECK_MPI(mpierr, "FILE_CLOSE!")
1474 : 0 : call xmpi_barrier(comm)
1475 : :
1476 [ # # ]: 0 : if (my_rank == 0) then
1477 : 0 : buf_size = 4
1478 : 0 : call wrtout(units, sjoin(" Reading norm of Coulomb vertex for testing purposes with ng:", itoa(buf_size)), newlines=1, pre_newlines=1)
1479 : 0 : call wrtout(units, " List of g-vectors in the Coulomb vertex:")
1480 [ # # ]: 0 : do ig=1,buf_size
1481 : 0 : call wrtout(units, sjoin(itoa(ig), ":", ltoa(m_gvec(:,ig))))
1482 : : end do
1483 [ # # # # : 0 : ABI_MALLOC(cwork, (buf_size))
# # # # #
# # # ]
1484 : 0 : call MPI_FILE_OPEN(xmpi_comm_self, cvx_filepath, MPI_MODE_RDONLY, xmpio_info, fh, mpierr)
1485 : 0 : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1486 : 0 : ierr = 0
1487 [ # # ]: 0 : band1_loop: do band1=1, ugb%nband_k
1488 [ # # ]: 0 : do band2=1, ugb%nband_k
1489 [ # # ]: 0 : ierr = ierr + 1; if (ierr == 6) exit band1_loop
1490 : 0 : offset = ((band2-1) * m_npw + (band1-1) * m_npw * ugb%nband_k) * xmpi_bsize_dpc
1491 : 0 : call MPI_FILE_READ_AT(fh, offset, cwork, buf_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1492 [ # # ]: 0 : ABI_HANDLE_MPIERR(mpierr)
1493 : 0 : call wrtout(units, sjoin(" For band1:", itoa(band1), ", band2:", itoa(band2)))
1494 [ # # # # ]: 0 : where (abs(cwork) < tol8)
1495 : : cwork = zero
1496 : : end where
1497 [ # # ]: 0 : write(msg, "(*(1x, es12.5))")abs(cwork(1:buf_size))
1498 : 0 : call wrtout(units, msg)
1499 : : end do
1500 : : end do band1_loop
1501 [ # # ]: 0 : ABI_FREE(cwork)
1502 : : end if
1503 : : #endif
1504 : :
1505 : : ! =============
1506 : : ! DEBUG SECTION
1507 : : ! =============
1508 : : #ifdef HAVE_MPI_IO
1509 : : if (my_rank == 0 .and. trust_no_one) then
1510 : : close(test_unt)
1511 : : if (open_file("test_mg", msg, newunit=test_unt, form="formatted", status="old", action="read") /= 0) then
1512 : : ABI_ERROR(msg)
1513 : : end if
1514 : :
1515 : : call MPI_FILE_OPEN(xmpi_comm_self, cvx_filepath, MPI_MODE_RDONLY, xmpio_info, fh, mpierr)
1516 : : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1517 : :
1518 : : ABI_MALLOC(cwork, (m_npw))
1519 : : max_abs_err = zero
1520 : : do ig=1, ugb%nband_k**2
1521 : : read(test_unt,*) band1, band2, ug12_batch(1:M_,1)
1522 : : offset = ((band2-1) * m_npw + (band1-1) * m_npw * ugb%nband_k) * xmpi_bsize_dpc
1523 : : call MPI_FILE_READ_AT(fh, offset, cwork, m_npw, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1524 : : ABI_HANDLE_MPIERR(mpierr)
1525 : :
1526 : : abs_err = maxval(abs(ug12_batch(1:M_,1) - cwork(1:M_)))
1527 : : max_abs_err = max(max_abs_err, abs_err)
1528 : : if (abs_err > zero) write(std_out, *)" For ig:", ig, "/", ugb%nband_k**2, "abs_err", abs_err
1529 : : !write(std_out, *)"1:", ug12_batch(1:M_,1); write(std_out, *)"2:", cwork(1:M_)
1530 : : end do
1531 : :
1532 : : close(test_unt)
1533 : : call MPI_FILE_CLOSE(fh, mpierr)
1534 : : ABI_CHECK_MPI(mpierr, "FILE_CLOSE!")
1535 : : ABI_FREE(cwork)
1536 : :
1537 : : write(std_out,*)" max_abs_err:", max_abs_err
1538 : : ABI_CHECK(max_abs_err < tol16, sjoin("max_abs_err:", ftoa(max_abs_err)))
1539 : : call wrtout(std_out, " Debugging section OK!!!")
1540 : : end if
1541 : : #endif
1542 : :
1543 : : ! Free memory
1544 : 0 : call uplan_1%free(); call uplan_2%free(); call uplan_m%free()
1545 : :
1546 [ # # ]: 0 : ABI_FREE(m_gvec)
1547 [ # # ]: 0 : ABI_FREE(ur1_batch)
1548 [ # # ]: 0 : ABI_FREE(ur2_batch)
1549 [ # # ]: 0 : ABI_FREE(ur12_batch)
1550 [ # # ]: 0 : ABI_FREE(ug12_batch)
1551 : 0 : ABI_FREE(sqrt_vc)
1552 : :
1553 [ # # ]: 0 : if (psps%usepaw == 1) then
1554 : 0 : call pawpwij_free(pwij)
1555 [ # # # # : 0 : ABI_FREE(pwij)
# # # # ]
1556 : 0 : call pawcprj_free(cprj1)
1557 [ # # # # : 0 : ABI_FREE(cprj1)
# # # # #
# ]
1558 [ # # ]: 0 : ABI_FREE(paw_rhotwg)
1559 : : end if
1560 : :
1561 : 0 : call cwtime_report(" cc4s_gamma", cpu, wall, gflops)
1562 : :
1563 [ # # # # : 0 : end subroutine cc4s_gamma
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1564 : : !!***
1565 : :
1566 : 0 : end module m_gwr_driver
1567 : : !!***
|
