Branch data Line data Source code
1 : : !!****m* ABINIT/m_wkk
2 : : !! NAME
3 : : !! m_wkk
4 : : !!
5 : : !! FUNCTION
6 : : !! Computation of the matrix elements of the screened interaction W_kk'.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt.
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_wkk
23 : :
24 : : !use, intrinsic :: iso_c_binding
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_xmpi
28 : : use m_mpinfo
29 : : use m_errors
30 : : use m_fft
31 : : use m_wfd
32 : : use m_hdr
33 : : use netcdf
34 : : use m_nctk
35 : : use m_dtset
36 : : use m_dtfil
37 : : use m_fstab
38 : : use m_ephtk
39 : :
40 : : use defs_abitypes, only : mpi_type
41 : : use defs_datatypes, only : pseudopotential_type
42 : : use m_gwdefs, only : GW_Q0_DEFAULT, cone_gw, czero_gw
43 : : use m_hide_blas, only : xgemm, xdotc
44 : : use m_special_funcs, only : gaussian
45 : : use m_time, only : cwtime, cwtime_report, timab, sec2str
46 : : use m_fstrings, only : itoa, ftoa, sjoin, ktoa, ltoa, strcat
47 : : use m_numeric_tools, only : arth, linspace ! print_arr
48 : : use m_io_tools, only : iomode_from_fname
49 : : use m_fftcore, only : ngfft_seq, sphereboundary, get_kg, print_ngfft
50 : : use m_crystal, only : crystal_t
51 : : use m_bz_mesh, only : kmesh_t, findqg0
52 : : use m_htetra, only : htetra_t
53 : : use m_gsphere, only : gsphere_t
54 : : use m_pawtab, only : pawtab_type
55 : : use m_io_screening, only : hscr_t, get_hscr_qmesh_gsph
56 : : use m_vcoul, only : vcoul_t
57 : : !use m_occ, only : get_fact_spin_tol_empty
58 : : use m_ebands, only : ebands_t, edos_t
59 : : use m_pstat, only : pstat_proc
60 : : use m_sigtk, only : sigtk_multiply_by_vc_sqrt
61 : : use m_screening, only : epsm1_t
62 : : use m_ddk, only : ddkop_t
63 : :
64 : : implicit none
65 : :
66 : : private
67 : : !!***
68 : :
69 : : public :: wkk_run ! Main entry point to compute Wkk' matrix elements.
70 : :
71 : : !----------------------------------------------------------------------
72 : :
73 : : contains !=====================================================
74 : : !!***
75 : :
76 : : !----------------------------------------------------------------------
77 : :
78 : : !!****f* m_wkk/wkk_run
79 : : !! NAME
80 : : !! wkk_run
81 : : !!
82 : : !! FUNCTION
83 : : !! Computation of the matrix elements of the screened interaction W_kk'
84 : : !! between two Cooper pairs.
85 : : !!
86 : : !! INPUTS
87 : : !! wfk0_path=String with the path to the GS unperturbed WFK file.
88 : : !! dtfil<datafiles_type>=Variables related to files.
89 : : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
90 : : !! dtset<dataset_type>=All input variables for this dataset.
91 : : !! cryst: Crystalline structure
92 : : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
93 : : !! wfk_hdr=Header of the WFK file.
94 : : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
95 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
96 : : !! comm=MPI communicator.
97 : : !!
98 : : !! OUTPUT
99 : : !! Results are written to file in netcdf format.
100 : : !!
101 : : !! SOURCE
102 : :
103 : 0 : subroutine wkk_run(wfk0_path, dtfil, ngfft, ngfftf, dtset, cryst, ebands, wfk_hdr, &
104 : 0 : pawtab, psps, mpi_enreg, comm)
105 : :
106 : : !Arguments ------------------------------------
107 : : !scalars
108 : : character(len=*),intent(in) :: wfk0_path
109 : : type(datafiles_type),intent(in) :: dtfil
110 : : type(dataset_type),intent(inout) :: dtset
111 : : type(crystal_t),intent(in) :: cryst
112 : : type(ebands_t),target,intent(in) :: ebands
113 : : type(hdr_type),intent(in) :: wfk_hdr
114 : : type(pseudopotential_type),intent(in) :: psps
115 : : type(mpi_type),intent(inout) :: mpi_enreg
116 : : integer,intent(in) :: comm
117 : : !arrays
118 : : integer,intent(in) :: ngfft(18), ngfftf(18)
119 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
120 : :
121 : : !Local variables ------------------------------
122 : : !scalars
123 : : integer,parameter :: LOG_MODQ = 1, LOG_MODK = 4, istwfk1 = 1, master = 0, ndat1 = 1
124 : : integer :: id_required, approx_type, ikxc, option_test, nkxc, ig
125 : : integer :: nkibz, my_rank, nsppol, iq_ibz, iq_bz, isym_qq, itim_qq
126 : : integer :: nspinor,nprocs, ii, spin, npw_x, npw_c ! max_npw_xc, min_npw_xc, my_is, cplex, ib,
127 : : integer :: bstart_k, bstop_k, nband_k, ncols, mpw, ierr !, ncerr
128 : : integer :: bstart_kp, bstop_kp, nband_kp, ik_bz, bmin, bmax, max_nb
129 : : integer :: ikp_ibz, isym_kp, trev_kp, npw_kp, istwf_kp, npw_kp_ibz, istwf_kp_ibz
130 : : integer :: ik_ibz, isym_k, trev_k, npw_k, istwf_k, npw_k_ibz, istwf_k_ibz
131 : : integer :: m_k, im_k, n_kp, in_kp, root_ncid, wkk_mode, ne
132 : : integer :: nfft,nfftf,mgfft,mgfftf,nkpg_kp,nkpg_k,edos_intmeth, nqlwl, scr_iomode ! cnt,
133 : : integer :: ikp_bz, my_ikp, my_nkp !, my_iq,
134 : : real(dp) :: cpu_all, wall_all, gflops_all, cpu_kp, wall_kp, gflops_kp ! cpu, wall, gflops,
135 : : real(dp) :: edos_step, edos_broad, e_mk, e_nkp, e_min, e_max, e_step, smear_mk, smear_nkp, faq
136 : : logical :: isirr_k, isirr_kp, qq_is_gamma, remove_exchange, print_time_kp
137 : 0 : type(wfd_t) :: wfd
138 [ # # # # : 0 : type(kmesh_t) :: qmesh, kmesh
# # # # ]
139 : 0 : type(gsphere_t),target :: gsph_x, gsph_c
140 [ # # ]: 0 : type(hscr_t),target :: hscr
141 [ # # # # : 0 : type(vcoul_t) :: vcp
# # # # #
# ]
142 : 0 : type(edos_t) :: edos
143 [ # # ]: 0 : type(epsm1_t) :: epsm1
144 [ # # # # : 0 : type(ddkop_t) :: ddkop
# # ]
145 : 0 : type(htetra_t) :: tetra
146 : : character(len=fnlen) :: screen_filepath
147 : : character(len=5000) :: msg !, qkp_string
148 : : !arrays
149 : : integer :: g0_k(3), g0_kp(3), g0_qq(3), units(2), work_ngfft(18), gmax(3), mapl_k(6), mapl_kp(6), mg0(3)
150 : 0 : integer,allocatable :: kg_kp(:,:), kg_k(:,:), gbound_kp(:,:), gbound_k(:,:), gbound_c(:,:), gbound_x(:,:)
151 : 0 : integer,allocatable :: nband(:,:), wfd_istwfk(:)
152 : 0 : integer, contiguous, pointer :: kg_c(:,:), kg_x(:,:)
153 : : real(dp) :: kk_ibz(3), kk_bz(3), kp_ibz(3), kp_bz(3), qq_bz(3), kk_diff(3) ! qq_ibz(3)
154 : : complex(gwp) :: ctmp_gwpc
155 : : character(len=fnlen) :: path
156 : : !arrays
157 : 0 : real(dp) :: n0(ebands%nsppol)
158 : 0 : real(dp),allocatable :: qlwl(:,:), kpg_kp(:,:), kpg_k(:,:), ug_kp(:,:,:), ug_k(:,:) !, cg_work(:,:)
159 : 0 : real(dp),allocatable :: work(:,:,:,:), e_mesh(:), e_args(:), wgt_mk(:), wgt_nkp(:,:) ! vcart_ibz(:,:,:,:),
160 : 0 : complex(gwp),allocatable :: cwork_ur(:), rhotwg_mn_x(:,:,:), rhotwg_mn_c(:,:,:), w_rhotwg_mn_c(:,:,:)
161 : 0 : complex(gwp),allocatable :: vc_sqrt_gx(:), ur_nkp(:,:), ur_mk(:,:), kxcg(:,:), mu_mn(:,:)
162 : 0 : complex(dp),allocatable :: w_ee(:,:) ! w_kkp(:,:,:,:)
163 : 0 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
164 : 0 : type(fstab_t),target,allocatable :: fstab(:)
165 : : !************************************************************************
166 : :
167 [ # # ]: 0 : ABI_CHECK(psps%usepaw == 0, "PAW not implemented")
168 [ # # ]: 0 : ABI_CHECK(dtset%nsppol == 1, "nsppol 2 implemented!")
169 [ # # ]: 0 : ABI_CHECK(dtset%nspinor == 1, "nspinor 2 implemented!")
170 : :
171 : : ABI_UNUSED(mpi_enreg%nproc)
172 : :
173 [ # # ]: 0 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm); units = [std_out, ab_out]
174 : 0 : call cwtime(cpu_all, wall_all, gflops_all, "start")
175 : :
176 : 0 : nsppol = ebands%nsppol; nspinor = ebands%nspinor
177 : :
178 : : ! Compute electron DOS.
179 [ # # ]: 0 : edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos
180 : : edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear
181 : 0 : edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha
182 : 0 : edos = ebands%get_edos(cryst, edos_intmeth, edos_step, edos_broad, comm)
183 : :
184 : : ! Get DOS per spin channel
185 [ # # ]: 0 : n0(:) = edos%gef(1:edos%nsppol)
186 [ # # ]: 0 : if (my_rank == master) then
187 : 0 : call edos%print(units)
188 [ # # ]: 0 : path = strcat(dtfil%filnam_ds(4), "_EDOS")
189 : 0 : call edos%write(path)
190 : : end if
191 : :
192 : : ! Find Fermi surface k-points.
193 : : ! TODO: support kptopt, change setup of k-points if tetra: fist tetra weights then k-points on the Fermi surface!
194 [ # # # # : 0 : ABI_MALLOC(fstab, (nsppol))
# # # # #
# # # #
# ]
195 : 0 : call fstab_init(fstab, ebands, cryst, dtset, tetra, comm)
196 : 0 : call tetra%free()
197 : :
198 : 0 : bmin = huge(1); bmax = -1
199 [ # # ]: 0 : do spin=1,nsppol
200 : 0 : bmin = min(bmin, fstab(spin)%bmin)
201 : 0 : bmax = max(bmax, fstab(spin)%bmax)
202 : : end do
203 : 0 : max_nb = bmax - bmin + 1
204 : :
205 : : ! Build linear mesh for W(e,e').
206 : 0 : spin = 1
207 : : associate (fs => fstab(spin))
208 [ # # # # : 0 : e_min = minval(ebands%eig(fs%bmin,:,spin)) - 0.1_dp * eV_Ha
# # # # ]
209 [ # # # # : 0 : e_max = maxval(ebands%eig(fs%bmax,:,spin)) + 0.1_dp * eV_Ha
# # # # ]
210 : 0 : e_step = 0.002 * eV_Ha
211 : 0 : ne = 1 + (e_max - e_min) / e_step
212 [ # # # # : 0 : ABI_CALLOC(w_ee, (ne, ne))
# # # # #
# # # # #
# # # # #
# ]
213 [ # # # # : 0 : ABI_MALLOC(e_mesh, (ne))
# # # # #
# # # ]
214 [ # # # # : 0 : ABI_MALLOC(e_args, (ne))
# # # # ]
215 [ # # # # : 0 : ABI_MALLOC(wgt_mk, (ne))
# # # # ]
216 [ # # # # : 0 : ABI_MALLOC(wgt_nkp, (ne, max_nb))
# # # # #
# # # # #
# # ]
217 [ # # # # : 0 : e_mesh = linspace(e_min, e_max, ne)
# # ]
218 : : end associate
219 : :
220 : : ! ================
221 : : ! HANDLE SCREENING
222 : : ! ================
223 : 0 : screen_filepath = dtfil%fnameabi_scr
224 [ # # ]: 0 : ABI_CHECK(dtfil%fnameabi_scr /= ABI_NOFILE, "SCR file must be specified in input!")
225 : :
226 : : ! Read g-sphere for correlation and qmesh from SCR file.
227 : 0 : call get_hscr_qmesh_gsph(screen_filepath, dtset, cryst, hscr, qmesh, gsph_c, qlwl, comm)
228 : 0 : call hscr%print(units, dtset%prtvol, header="Header of the SCR file")
229 : :
230 : 0 : nqlwl = size(qlwl, dim=2)
231 [ # # ]: 0 : if (nqlwl == 0) then
232 : 0 : nqlwl=1
233 [ # # # # ]: 0 : ABI_MALLOC(qlwl,(3,nqlwl))
234 [ # # ]: 0 : qlwl(:,nqlwl)= GW_Q0_DEFAULT
235 : : write(msg,'(3a,i0,a,3f9.6)')&
236 : 0 : "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,&
237 : 0 : "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1)
238 : 0 : ABI_COMMENT(msg)
239 : : end if
240 : :
241 : : ! Init g-sphere for the exchange part from ecutsigx.
242 : 0 : call gsph_c%extend(cryst, dtset%ecutsigx, gsph_x)
243 : :
244 : : ! Initialize Coulomb term on the IBZ of the qmesh. Use largest G-sphere.
245 : 0 : call kmesh%init(cryst, wfk_hdr%nkpt, wfk_hdr%kptns, dtset%kptopt)
246 : 0 : faq = one / (cryst%ucvol*kmesh%nbz)
247 : :
248 : : ! TODO:
249 : : ! Here we sort the pp_mesh by stars so that we can split the pp wavevectors in blocks and therefore
250 : : ! reduce the number of wavevectors in the IBZ that must be stored in memory.
251 : : !call qmesh%pack_by_stars()
252 : :
253 : 0 : npw_x = gsph_x%ng; npw_c = gsph_c%ng !; min_npw_xc = min(npw_x, npw_c); max_npw_xc = max(npw_x, npw_c)
254 [ # # ]: 0 : if (gsph_x%ng >= gsph_c%ng) then
255 : : call vcp%init(gsph_x, cryst, qmesh, kmesh, dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_x%ng, &
256 : 0 : nqlwl, qlwl, comm)
257 : : else
258 : : call vcp%init(gsph_c, cryst, qmesh, kmesh, dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_c%ng, &
259 : 0 : nqlwl, qlwl, comm)
260 : : end if
261 : :
262 [ # # ]: 0 : ABI_FREE(qlwl)
263 : :
264 [ # # ]: 0 : if (my_rank == master) then
265 : 0 : call kmesh%print(units, header="K-mesh for wavefunctions", prtvol=dtset%prtvol)
266 : 0 : call gsph_x%print(units, dtset%prtvol, header="G-sphere for exchange")
267 : 0 : call gsph_c%print(units, dtset%prtvol, header="G-sphere for correlation")
268 : 0 : call vcp%print(units, prtvol=dtset%prtvol)
269 : : end if
270 : 0 : call kmesh%free()
271 : :
272 [ # # ]: 0 : ABI_CHECK_IGE(npw_x, 1, "npw_x <= 1")
273 [ # # ]: 0 : ABI_CHECK_IGE(npw_c, 1, "npw_c <= 1")
274 : :
275 : : ! Note that in this case, the sphere is always Gamma-centered i.e. it does not depend on the qq wavevector
276 : 0 : kg_c => gsph_c%gvec(:, 1:npw_c)
277 : 0 : kg_x => gsph_x%gvec(:, 1:npw_x)
278 : :
279 : : ! Initialize the wave function descriptor.
280 : : ! Only wavefunctions on the FS are stored in wfd.
281 : : ! Need all k-points on the FS because of k+q, spin is not distributed for the time being.
282 : : ! It would be possible to reduce the memory allocated per MPI-rank via OpenMP.
283 : 0 : nkibz = wfk_hdr%nkpt
284 [ # # # # : 0 : ABI_MALLOC(nband, (nkibz, nsppol))
# # # # #
# # # # #
# # ]
285 [ # # # # : 0 : ABI_MALLOC(bks_mask, (dtset%mband, nkibz, nsppol))
# # # # #
# # # # #
# # # # #
# ]
286 [ # # # # : 0 : ABI_MALLOC(keep_ur, (dtset%mband, nkibz, nsppol))
# # # # #
# # # # #
# # ]
287 [ # # # # : 0 : nband = dtset%mband; bks_mask = .False.; keep_ur = .False.
# # # # #
# # # # #
# # ]
288 : :
289 [ # # ]: 0 : do spin=1,ebands%nsppol
290 : 0 : associate (fs => fstab(spin))
291 [ # # ]: 0 : do ik_bz=1,fs%nkfs
292 : 0 : ik_ibz = fs%indkk_fs(1, ik_bz)
293 : 0 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz)
294 [ # # ]: 0 : bks_mask(bstart_k:bstart_k+nband_k-1, ik_ibz, spin) = .True.
295 : : end do
296 : : end associate
297 : : end do
298 : :
299 : : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
300 : : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
301 [ # # # # : 0 : ABI_MALLOC(wfd_istwfk, (nkibz))
# # # # #
# # # ]
302 [ # # ]: 0 : wfd_istwfk = 1
303 : :
304 : : call wfd%init(cryst, pawtab, psps, keep_ur, dtset%mband, nband, nkibz, nsppol, bks_mask,&
305 : : dtset%nspden, nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,&
306 : 0 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
307 : :
308 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
309 : :
310 [ # # ]: 0 : ABI_FREE(nband)
311 [ # # ]: 0 : ABI_FREE(bks_mask)
312 [ # # ]: 0 : ABI_FREE(keep_ur)
313 [ # # ]: 0 : ABI_FREE(wfd_istwfk)
314 : :
315 : : ! Read wavefunctions.
316 : 0 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
317 : :
318 : : ! TODO
319 : : ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
320 : : ! We also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
321 : : ! that will be used to symmetrize the wavefunctions in G-space.
322 : : ! This is the maximum number of PWs for all possible k+q treated.
323 : :
324 : 0 : call ephtk_get_mpw_gmax(nkibz, wfk_hdr%kptns, dtset%ecut, cryst%gmet, mpw, gmax, comm)
325 : : ! FIXME
326 : : !gmax = 2 * gmax
327 : :
328 : 0 : call ddkop%init(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft)
329 : :
330 : : #if 0
331 : : call cwtime(cpu, wall, gflops, "start", msg=" Computing v_nk matrix elements for all states on the FS...")
332 : :
333 : : ABI_CALLOC(vcart_ibz, (3, bmin:bmax, nkibz, nsppol))
334 : : ABI_MALLOC(cg_work, (2, mpw * nspinor))
335 : :
336 : : cnt = 0
337 : : do spin=1,nsppol
338 : : associate (fs => fstab(spin))
339 : : do ik_ibz=1,ebands%nkpt
340 : : kk_ibz = ebands%kptns(:, ik_ibz)
341 : : npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz)
342 : : ! NB: The two checks below are global --> all procs will cycle.
343 : : if (all(bks_mask(:, ik_ibz, spin) .eqv. .False.)) cycle
344 : : if (npw_k == 1) cycle
345 : : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
346 : :
347 : : call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk_ibz, istwf_k, npw_k, wfd%kdata(ik_ibz)%kg_k)
348 : :
349 : : do band_k=fs%bmin,fs%bmax
350 : : if (.not. bks_mask(band_k, ik_ibz, spin)) cycle
351 : : call wfd%copy_cg(band_k, ik_ibz, spin, cg_work)
352 : : eig0nk = ebands%eig(band_k, ik_ibz, spin)
353 : : vk = ddkop%get_vdiag(eig0nk, istwf_k, npw_k, wfd%nspinor, cg_work, cwaveprj0)
354 : : vcart_ibz(:, band_k, ik_ibz, spin) = vk
355 : : end do
356 : : end do
357 : : end associate
358 : : end do ! spin
359 : :
360 : : call xmpi_sum(vcart_ibz, comm, ierr)
361 : : ABI_FREE(cg_work)
362 : : ABI_FREE(vcart_ibz)
363 : : call cwtime_report(" Velocities", cpu, wall, gflops)
364 : : #endif
365 : 0 : call ddkop%free()
366 : :
367 : 0 : call ngfft_seq(work_ngfft, gmax)
368 [ # # # # : 0 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
# # # # #
# # # # #
# # # # #
# ]
369 : :
370 : : ! FFT meshes from input file, not necessarily equal to the ones found in the external files.
371 : : ! NB: ur arrays are always allocated with nfft and not with product(ngfft(4:6)).
372 [ # # ]: 0 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
373 [ # # # # ]: 0 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
374 : :
375 [ # # ]: 0 : call print_ngfft([std_out], ngfft, header="FFT mesh")
376 : :
377 : : ! Set the FFT mesh
378 : 0 : call wfd%change_ngfft(cryst, psps, ngfft)
379 : 0 : call wfd%print(units, header="Wavefunctions for GWPT calculation.")
380 : :
381 [ # # # # : 0 : ABI_MALLOC(gbound_kp, (2*mgfft+8, 2))
# # # # #
# # # # #
# # ]
382 [ # # # # : 0 : ABI_MALLOC(gbound_k, (2*mgfft+8, 2))
# # # # ]
383 [ # # # # : 0 : ABI_MALLOC(gbound_c, (2*mgfft+8, 2))
# # # # ]
384 [ # # # # : 0 : ABI_MALLOC(gbound_x, (2*mgfft+8, 2))
# # # # ]
385 : :
386 : : ! ====================================
387 : : ! This is the g-sphere for W_{gg'}(qq)
388 : : ! ====================================
389 : : ! Init g-sphere for the exchange part from ecutsigx.
390 : 0 : call sphereboundary(gbound_c, istwfk1, kg_c, mgfft, npw_c)
391 : 0 : call sphereboundary(gbound_x, istwfk1, kg_x, mgfft, npw_x)
392 : :
393 [ # # # # : 0 : ABI_MALLOC(vc_sqrt_gx, (npw_x))
# # # # #
# # # ]
394 : :
395 : : ! Allocate g-vectors centered on k, k'
396 [ # # # # : 0 : ABI_MALLOC(kg_k, (3, mpw))
# # # # #
# # # ]
397 [ # # # # : 0 : ABI_MALLOC(kg_kp, (3, mpw))
# # # # ]
398 : :
399 : 0 : spin = 1
400 : : associate (fs => fstab(spin))
401 [ # # # # : 0 : ABI_MALLOC_OR_DIE(ur_nkp, (nfft*nspinor, fs%bmax-fs%bmin+1), ierr)
# # # # #
# # # # #
# # # # ]
402 [ # # # # : 0 : ABI_MALLOC_OR_DIE(ur_mk, (nfft*nspinor, fs%bmax-fs%bmin+1), ierr)
# # # # #
# # # # #
# # # # #
# ]
403 [ # # # # : 0 : ABI_MALLOC(cwork_ur, (nfft*nspinor))
# # # # #
# # # ]
404 : :
405 : : ! Read symmetrized em1 from file
406 : 0 : id_required = 4; ikxc = 0; approx_type = 0; option_test = 0; nkxc = 0
407 [ # # # # : 0 : ABI_MALLOC(kxcg, (nfftf, nkxc))
# # # # ]
408 : 0 : scr_iomode = iomode_from_fname(screen_filepath)
409 : :
410 : 0 : call epsm1%from_file(screen_filepath, qmesh%nibz, npw_c, comm)
411 : : call epsm1%mkdump(vcp, npw_c, kg_c, nkxc, kxcg, id_required, approx_type, &
412 : 0 : ikxc, option_test, dtfil%fnameabo_scr, scr_iomode, nfftf, ngfft, comm)
413 : :
414 : 0 : call epsm1%malloc_epsm1_qbz(npw_c, epsm1%nomega)
415 [ # # ]: 0 : ABI_FREE(kxcg)
416 : :
417 : 0 : call pstat_proc%print(_PSTAT_ARGS_)
418 : :
419 : 0 : wkk_mode = 1
420 : : ! Master creates the netcdf file used to store the results of the calculation.
421 [ # # ]: 0 : if (my_rank == master) then
422 [ # # ]: 0 : NCF_CHECK(nctk_open_create(root_ncid, strcat(dtfil%filnam_ds(4), "_WKK.nc") , xmpi_comm_self))
423 [ # # ]: 0 : NCF_CHECK(cryst%ncwrite(root_ncid))
424 [ # # ]: 0 : NCF_CHECK(ebands%ncwrite(root_ncid))
425 [ # # ]: 0 : NCF_CHECK(edos%ncwrite(root_ncid))
426 : :
427 : : !ncerr = nctk_def_dims(ncid, [ &
428 : : ! nctkdim_t("smat_bsize1", smat_bsize1), nctkdim_t("smat_bsize2", smat_bsize2) &
429 : : ! ], defmode=.True.)
430 : : !NCF_CHECK(ncerr)
431 : : !ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
432 : : ! "gwr_completed", "scf_iteration" &
433 : : !])
434 : : !NCF_CHECK(ncerr)
435 : : !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
436 : : ! "wr_step", "gwr_boxcutmin", "cosft_duality_error", "regterm" &
437 : : !])
438 : : !NCF_CHECK(ncerr)
439 : :
440 : : !ncerr = nctk_def_arrays(ncid, [ &
441 : : ! nctkarr_t("gwr_task", "char", "character_string_length"), &
442 : : !])
443 : : !NCF_CHECK(ncerr)
444 : : !NCF_CHECK(nctk_set_datamode(ncid))
445 : :
446 : : !ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
447 : : ! "gwr_completed", "sig_diago", "b1gw", "b2gw", "symsigma", "symchi", "scf_iteration"], &
448 : : ! [0, merge(1, 0, gwr%sig_diago), gwr%b1gw, gwr%b2gw, gwr%dtset%symsigma, dtset%symchi, gwr%scf_iteration])
449 : : !NCF_CHECK(ncerr)
450 : :
451 : : !ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
452 : : ! "wr_step", "ecuteps", "ecut", "ecutwfn", "ecutsigx", "gwr_boxcutmin", &
453 : : ! "cosft_duality_error", "regterm"], &
454 : : ! [gwr%wr_step, dtset%ecuteps, dtset%ecut, dtset%ecutwfn, dtset%ecutsigx, dtset%gwr_boxcutmin, &
455 : : ! gwr%ft_max_error(1), gwr%ft_max_error(2), gwr%ft_max_error(3), gwr%cosft_duality_error, regterm &
456 : : ! ])
457 : : !NCF_CHECK(ncerr)
458 : : !NCF_CHECK(nf90_put_var(ncid, vid("gwr_task"), trim(dtset%gwr_task)))
459 : :
460 [ # # ]: 0 : NCF_CHECK(nf90_close(root_ncid))
461 : : end if
462 : 0 : call xmpi_barrier(comm)
463 : :
464 : : ! Open WKK.nc file and go to data mode.
465 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(root_ncid, strcat(dtfil%filnam_ds(4), "_WKK.nc") , comm))
466 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(root_ncid))
467 : :
468 : 0 : mg0 = [1, 1, 1]
469 : 0 : remove_exchange = .True.
470 : :
471 [ # # # # : 0 : ABI_CALLOC(mu_mn, (bmin:bmax, bmin:bmax))
# # # # #
# # # # #
# # # # #
# ]
472 : :
473 : : ! Loop over k'-points in the energy window.
474 [ # # # # ]: 0 : do ikp_bz=1,fs%nkfs
475 : 0 : my_ikp = ikp_bz; my_nkp = fs%nkfs
476 : :
477 [ # # # # : 0 : print_time_kp = my_rank == 0 .and. (ikp_bz <= LOG_MODK .or. mod(ikp_bz, LOG_MODK) == 0)
# # ]
478 : 0 : if (print_time_kp) call cwtime(cpu_kp, wall_kp, gflops_kp, "start")
479 : :
480 : : ! The k-point and the symmetries relating the BZ k-point to the IBZ.
481 [ # # ]: 0 : kp_bz = fs%kpts(:, ikp_bz)
482 : 0 : ikp_ibz = fs%indkk_fs(1, ikp_bz) ; isym_kp = fs%indkk_fs(2, ikp_bz)
483 [ # # ]: 0 : trev_kp = fs%indkk_fs(6, ikp_bz); g0_kp = fs%indkk_fs(3:5,ikp_bz)
484 [ # # # # : 0 : isirr_kp = (isym_kp == 1 .and. trev_kp == 0 .and. all(g0_kp == 0))
# # # # ]
485 [ # # ]: 0 : mapl_kp = fs%indkk_fs(:, ikp_bz)
486 : :
487 [ # # ]: 0 : kp_ibz = ebands%kptns(:,ikp_ibz)
488 : 0 : istwf_kp_ibz = wfd%istwfk(ikp_ibz); npw_kp_ibz = wfd%npwarr(ikp_ibz)
489 : :
490 : : ! Get npw_kp, kg_kp for this kp_bz.
491 : : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kp_bz, ikp_ibz, isirr_kp, dtset%nloalg, & ! in
492 [ # # ]: 0 : istwf_kp, npw_kp, kg_kp, nkpg_kp, kpg_kp, gbound_kp) ! out
493 : :
494 : : ! Rotate from IBZ to BZ and compute ur_nkp for all n_kp bands.
495 : 0 : bstart_kp = fs%bstart_cnt_ibz(1, ikp_ibz); nband_kp = fs%bstart_cnt_ibz(2, ikp_ibz); bstop_kp = bstart_kp + nband_kp - 1
496 [ # # # # : 0 : ABI_CALLOC(ug_kp, (2, npw_kp*nspinor, nband_kp))
# # # # #
# # # # #
# # # # #
# # # ]
497 : : call wfd%rotate_cg(bstart_kp, nband_kp, spin, kp_ibz, npw_kp, kg_kp, istwf_kp, &
498 : 0 : cryst, mapl_kp, gbound_kp, work_ngfft, work, ug_kp, urs_kbz=ur_nkp)
499 : 0 : do ii=1, nband_kp
500 : : !print *, "isirr_kp:", isirr_kp
501 : : !print *, "g:", sum(ug_kp(1,:,ii)**2 + ug_kp(2,:,ii)**2)
502 : : !print *, "r:", sum(abs(ur_nkp(:,ii)) ** 2) / nfft ! * cryst%ucvol /
503 : : end do
504 : :
505 : : ! Loop over k-points in the energy window.
506 [ # # ]: 0 : do ik_bz=1,fs%nkfs
507 : :
508 [ # # ]: 0 : mapl_k = fs%indkk_fs(:, ik_bz)
509 [ # # ]: 0 : ik_ibz = mapl_k(1); isym_k = mapl_k(2); trev_k = mapl_k(6); g0_k = mapl_k(3:5)
510 [ # # # # : 0 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
# # # # ]
511 [ # # ]: 0 : kk_ibz = ebands%kptns(:, ik_ibz)
512 : 0 : istwf_k_ibz = wfd%istwfk(ik_ibz); npw_k_ibz = wfd%npwarr(ik_ibz)
513 : :
514 [ # # ]: 0 : kk_bz = fs%kpts(:, ik_bz)
515 [ # # ]: 0 : kk_ibz = ebands%kptns(:,ik_ibz)
516 : :
517 : : ! Get npw_k, kg_k for this k.
518 : : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kk_bz, ik_ibz, isirr_k, dtset%nloalg, & ! in
519 [ # # ]: 0 : istwf_k, npw_k, kg_k, nkpg_k, kpg_k, gbound_k) ! out
520 [ # # ]: 0 : ABI_FREE(kpg_k)
521 : :
522 : 0 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz); bstop_k = bstart_k + nband_k - 1
523 [ # # # # : 0 : ABI_MALLOC(ug_k, (2, npw_k*nspinor))
# # # # #
# # # ]
524 [ # # # # : 0 : ABI_MALLOC(rhotwg_mn_x, (npw_x*nspinor, nband_k, nband_kp))
# # # # #
# # # # #
# # # # #
# ]
525 : :
526 [ # # ]: 0 : do m_k=bstart_k, bstop_k
527 : : ! Rotate from IBZ to BZ and compute ur_mk
528 : 0 : im_k = m_k - bstart_k + 1
529 : : call wfd%rotate_cg(m_k, ndat1, spin, kk_ibz, npw_k, kg_k, istwf_k, &
530 : 0 : cryst, mapl_k, gbound_k, work_ngfft, work, ug_k, urs_kbz=ur_mk(:,im_k))
531 : :
532 [ # # ]: 0 : do n_kp=bstart_kp, bstop_kp
533 : 0 : in_kp = n_kp - bstart_kp + 1
534 [ # # # # ]: 0 : cwork_ur = conjg(ur_nkp(:,in_kp)) * ur_mk(:,im_k)
535 : 0 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_mn_x(:, im_k, in_kp))
536 : 0 : call sigtk_multiply_by_vc_sqrt("N", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_mn_x(:, im_k, in_kp))
537 : : end do
538 : : end do
539 : :
540 [ # # ]: 0 : ABI_FREE(ug_k)
541 : :
542 : : ! Identify q and G0 where q + G0 = k_GW - k_i
543 [ # # ]: 0 : kk_diff = kk_bz - kp_bz
544 : :
545 : 0 : call findqg0(iq_bz, g0_qq, kk_diff, qmesh%nbz, qmesh%bz, mG0)
546 : :
547 : : ! Find the corresponding irred qq-point in qmesh.
548 : 0 : call qmesh%get_bz_item(iq_bz, qq_bz, iq_ibz, isym_qq, itim_qq)
549 : 0 : qq_is_gamma = sum(qq_bz**2) < tol14
550 : :
551 : : ! Find k + q in the extended zone and extract symmetry info.
552 : : ! Be careful here because there are two umklapp vectors to be considered as:
553 : : !
554 : : ! k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
555 : : !
556 : : ! Get Fourier components of the Coulomb interaction in the BZ
557 : : ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2
558 : : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
559 : : ! NOTE: vc_sqrt_gx is dimensioned with npw_x --> use rottb table from gsph_x.
560 [ # # ]: 0 : do ig=1,npw_x
561 : 0 : vc_sqrt_gx(gsph_x%rottb(ig, itim_qq, isym_qq)) = vcp%vc_sqrt(ig, iq_ibz)
562 : : end do
563 : :
564 : 0 : call epsm1%rotate_iqbz(iq_bz, hscr%nomega, npw_c, gsph_c, qmesh, remove_exchange)
565 : :
566 : : ! Compute \sum_{gg'} (M_g^{kn,k'm})^* W_{gg'}(q) M_g'^{kn,k'm} with k' = k + q.
567 : : ! TODO: nspinor
568 [ # # # # : 0 : ABI_MALLOC(rhotwg_mn_c, (npw_c*nspinor, nband_k, nband_kp))
# # # # #
# # # # #
# # # # #
# ]
569 [ # # # # : 0 : rhotwg_mn_c(:,:,:) = rhotwg_mn_x(1:npw_c,:,:)
# # ]
570 : :
571 : 0 : ncols = nband_k * nband_kp
572 [ # # # # : 0 : ABI_MALLOC(w_rhotwg_mn_c, (npw_c*nspinor, nband_k, nband_kp))
# # # # #
# # # # #
# # # # #
# ]
573 : :
574 : : call xgemm("N", "N", npw_c, ncols, npw_c, cone_gw, epsm1%epsm1_qbz(:,:,1), npw_c, rhotwg_mn_c(:,:,1), npw_c*nspinor, &
575 [ # # # # : 0 : czero_gw, w_rhotwg_mn_c(:,:,1), npw_c)
# # # # #
# ]
576 : :
577 : : !do in_kp=1, nband_k
578 : : ! do im_k=1, nband_kq
579 : : ! w_rhotwg_mn_c(:,in_kp,im_k) = matmul(epsm1%epsm1_qbz(:,:,1), rhotwg_mn_c(:,im_k,in_kp))
580 : : ! end do
581 : : !end do
582 : : !print *, " rhotwg_mn_c", maxval(abs(rhotwg_mn_c)), maxloc(abs(rhotwg_mn_c))
583 : : !print *, "w_rhotwg_mn_c:", maxval(abs(w_rhotwg_mn_c))
584 : :
585 : : ! Precompute delta(e - e_nkp - ef).
586 [ # # ]: 0 : do n_kp=bstart_kp, bstop_kp
587 : 0 : in_kp = n_kp - bstart_kp + 1
588 : 0 : e_nkp = ebands%eig(n_kp, ik_ibz, spin)
589 [ # # # # ]: 0 : e_args = e_mesh - e_nkp
590 : 0 : smear_nkp = 0.1_dp * eV_Ha
591 [ # # ]: 0 : wgt_nkp(:, im_k) = gaussian(e_args, smear_nkp)
592 : : end do
593 : :
594 [ # # ]: 0 : do n_kp=bstart_kp, bstop_kp
595 : 0 : in_kp = n_kp - bstart_kp + 1
596 : 0 : e_nkp = ebands%eig(n_kp, ikp_ibz, spin)
597 : : !fs%tetra_wtk(n_kp - fs%bmin + 1, ikp_ibz)
598 [ # # ]: 0 : do m_k=bstart_k, bstop_k
599 : 0 : im_k = m_k - bstart_k + 1
600 : 0 : e_mk = ebands%eig(m_k, ik_ibz, spin)
601 : : !fs%tetra_wtk(m_k - fs%bmin + 1, ik_ibz)
602 : 0 : smear_mk = 0.1_dp * eV_Ha
603 : 0 : ctmp_gwpc = xdotc(npw_c*nspinor, rhotwg_mn_c(:,im_k,in_kp), 1, w_rhotwg_mn_c(:,im_k,in_kp), 1)
604 : : if (remove_exchange) then
605 : 0 : ctmp_gwpc = ctmp_gwpc + xdotc(npw_x*nspinor, rhotwg_mn_x(:,im_k,in_kp), 1, rhotwg_mn_x(:,im_k,in_kp), 1)
606 : : end if
607 : 0 : ctmp_gwpc = faq * ctmp_gwpc
608 : : !print *, "ctmp_gwpc:", ctmp_gwpc
609 : : !mu_mn(m_k, n_kp) = mu_mn(m_k, n_kp) + &
610 : : ! ctmp_gwpc * gaussian(e_mk - ebands%fermie, smear_mk) * gaussian(e_mp, smear_mk)
611 : : ! Computes: A := A + alpha * x * y^T
612 : : !w_ee =
613 : : !call dger(ne, ne, alpha, x, incx, y, incy, w_ee, lda)
614 : : end do
615 : : end do
616 : :
617 [ # # ]: 0 : ABI_FREE(rhotwg_mn_x)
618 [ # # ]: 0 : ABI_FREE(rhotwg_mn_c)
619 [ # # ]: 0 : ABI_FREE(w_rhotwg_mn_c)
620 : : end do ! ik_bz
621 : :
622 [ # # ]: 0 : if (print_time_kp) then
623 : 0 : call inds2str(1, "My kp-point", my_ikp, my_nkp, fs%nkfs, msg)
624 [ # # ]: 0 : call cwtime_report(msg, cpu_kp, wall_kp, gflops_kp); if (ikp_bz == LOG_MODK) call wrtout(std_out, "...", do_flush=.True.)
625 : : end if
626 : :
627 [ # # ]: 0 : ABI_FREE(ug_kp)
628 [ # # ]: 0 : ABI_FREE(kpg_kp)
629 : : end do ! ikp_bz
630 : : end associate
631 : :
632 [ # # ]: 0 : ABI_FREE(ur_nkp)
633 [ # # ]: 0 : ABI_FREE(ur_mk)
634 : : !end do ! my_is
635 : :
636 : 0 : call cwtime_report(" wkk calculation", cpu_all, wall_all, gflops_all, end_str=ch10)
637 : :
638 : 0 : call xmpi_sum(w_ee, comm, ierr)
639 : 0 : call xmpi_sum(mu_mn, comm, ierr)
640 : :
641 : : if (my_rank == master) then
642 : : !do m_k=bmin, bmax
643 : : ! do n_kp=bmin, bmax
644 : : ! call print_arr(units, mu_mn, max_r=max_nb, max_c=max_nb)
645 : : ! end do
646 : : !end do
647 : : !call wrtout(units, sjoin("Re mu: ", ftoa(sum(real(mu_mn)))))
648 : : !call wrtout(units, sjoin("Im mu: ", ftoa(sum(aimag(mu_mn)))))
649 : : end if
650 : :
651 : : ! Set vqk_completed to 1 so that we can easily check if restarted is needed.
652 : : !if (wkk_mode /= 0) then
653 : : !if (my_rank == master) then
654 : : ! NCF_CHECK(nf90_put_var(root_ncid, root_vid("vqk_completed"), 1))
655 : : !end if
656 [ # # ]: 0 : NCF_CHECK(nf90_close(root_ncid))
657 : 0 : call xmpi_barrier(comm)
658 : : !end if
659 : :
660 : : ! Output some of the results to ab_out for testing purposes
661 : : !call gstore%print_for_abitests(dtset)
662 : :
663 : : ! Free memory
664 [ # # ]: 0 : ABI_FREE(kg_kp)
665 [ # # ]: 0 : ABI_FREE(kg_k)
666 [ # # ]: 0 : ABI_FREE(cwork_ur)
667 [ # # ]: 0 : ABI_FREE(work)
668 [ # # ]: 0 : ABI_FREE(gbound_kp)
669 [ # # ]: 0 : ABI_FREE(gbound_k)
670 [ # # ]: 0 : ABI_FREE(gbound_c)
671 [ # # ]: 0 : ABI_FREE(gbound_x)
672 [ # # ]: 0 : ABI_FREE(vc_sqrt_gx)
673 [ # # ]: 0 : ABI_FREE(w_ee)
674 [ # # ]: 0 : ABI_FREE(e_mesh)
675 [ # # ]: 0 : ABI_FREE(e_args)
676 [ # # ]: 0 : ABI_FREE(wgt_mk)
677 [ # # ]: 0 : ABI_FREE(wgt_nkp)
678 [ # # ]: 0 : ABI_FREE(mu_mn)
679 : :
680 : 0 : call wfd%free(); call vcp%free(); call qmesh%free(); call gsph_x%free(); call gsph_c%free(); call hscr%free(); call edos%free()
681 : 0 : call epsm1%free()
682 : :
683 [ # # ]: 0 : do spin=1,ebands%nsppol
684 : 0 : call fstab(spin)%free()
685 : : end do
686 [ # # # # : 0 : ABI_FREE(fstab)
# # # # #
# # # # #
# # # # #
# # # #
# ]
687 : :
688 : 0 : call xmpi_barrier(comm) ! This to make sure that the parallel output of WKK is completed
689 [ # # # # : 0 : call cwtime_report(" wkk_run: MPI barrier before returning.", cpu_all, wall_all, gflops_all, end_str=ch10, comm=comm)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
690 : :
691 : : contains
692 : :
693 : 0 : subroutine inds2str(ilevel, prefix, my_ik, my_nk, nk_tot, out_str)
694 : : character(len=*),intent(in) :: prefix
695 : : integer,intent(in) :: ilevel, my_ik, my_nk, nk_tot
696 : : character(len=*),intent(out) :: out_str
697 : :
698 [ # # # # ]: 0 : out_str = sjoin(prefix, itoa(my_ik), "/", itoa(my_nk), "[", itoa(nk_tot), "]")
699 [ # # # # : 0 : out_str = repeat(' ', 4 * ilevel) // trim(out_str)
# # ]
700 : 0 : end subroutine inds2str
701 : :
702 : : integer function root_vid(var_name)
703 : : character(len=*),intent(in) :: var_name
704 : : root_vid = nctk_idname(root_ncid, var_name)
705 : : end function root_vid
706 : :
707 : : !integer function spin_vid(var_name)
708 : : ! character(len=*),intent(in) :: var_name
709 : : ! spin_vid = nctk_idname(spin_ncid, var_name)
710 : : !end function spin_vid
711 : :
712 : : end subroutine wkk_run
713 : : !!***
714 : :
715 : : end module m_wkk
716 : : !!***
|
