Branch data Line data Source code
1 : : !!****m* ABINIT/m_conducti
2 : : !! NAME
3 : : !! m_conducti
4 : : !!
5 : : !! FUNCTION
6 : : !! This program computes the elements of the optical frequency dependent
7 : : !! conductivity tensor and the conductivity along the three principal axes
8 : : !! from the Kubo-Greenwood formula.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2002-2026 ABINIT group (VRecoules, PGhosh, SMazevet, SM, SVinko, NBrouwer)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : ! nvtx related macro definition
25 : : #include "nvtx_macros.h"
26 : :
27 : : module m_conducti
28 : :
29 : : use defs_basis
30 : : use m_errors
31 : : use m_abicore
32 : : use m_xmpi
33 : : use m_wfk
34 : : use m_hdr
35 : : use m_nctk
36 : : use netcdf
37 : :
38 : : use defs_abitypes, only : MPI_type
39 : : use m_io_tools, only : open_file, close_unit, get_unit
40 : : use m_fstrings, only : sjoin
41 : : use m_matrix, only : matr3inv
42 : : use m_hide_lapack, only : jacobi
43 : : use m_occ, only : getnel
44 : : use m_geometry, only : metric
45 : : use m_splines, only : splint,spline
46 : : use m_mpinfo, only : distrb2,init_mpi_enreg,destroy_mpi_enreg,proc_distrb_cycle,initmpi_band
47 : : use m_time, only : timab
48 : : #if defined(HAVE_GPU_MARKERS)
49 : : use m_nvtx
50 : : use m_nvtx_data
51 : : #endif
52 : :
53 : : implicit none
54 : :
55 : : private
56 : : !!***
57 : :
58 : : public :: conducti_paw
59 : : public :: conducti_paw_core
60 : : public :: conducti_nc
61 : :
62 : : !I/O parameters
63 : : !Set to true to use netcdf-MPIIO when available
64 : : logical,parameter :: use_netcdf_mpiio=.true.
65 : :
66 : : !Number of (n,m) band pairs buffered before applying the Kubo-Greenwood formula
67 : : !(pair data: 10 reals and 1 logical per pair, about 2.7 MB, to stay in cache)
68 : : integer,parameter :: kg_pairs_blk=32768
69 : : !Same with GPU: larger buffer (about 88 MB) to limit the number of transfers and kernel launches
70 : : integer,parameter :: kg_pairs_blk_gpu=1048576
71 : : !With GPU, each frequency sums over the pairs by (at most) kg_gpu_nchunk chunks of pairs,
72 : : ! treated by different threads; the partial sums are added afterwards
73 : : integer,parameter :: kg_gpu_nchunk=128
74 : :
75 : : !Maximum memory used to store the dipoles of a block of bands read with MPI-IO
76 : : real(dp),parameter :: dipoles_blk_mem=512._dp*1024._dp**2
77 : :
78 : : !Number of bands buffered before computing the core-valence contributions
79 : : integer,parameter :: core_bands_blk=1024
80 : :
81 : : !!***
82 : :
83 : : contains
84 : : !!***
85 : :
86 : : !!****f* m_conducti/conducti_paw
87 : : !! NAME
88 : : !! conducti_paw
89 : : !!
90 : : !! FUNCTION
91 : : !! This program computes the elements of the optical frequency dependent
92 : : !! conductivity tensor and the conductivity along the three principal axes
93 : : !! from the Kubo-Greenwood formula for PAW formalism
94 : : !!
95 : : !! INPUTS
96 : : !! filnam=generic name for input data
97 : : !! filnam_out=generic name for output data
98 : : !! [varocc]=if true, read arbitrary occupations from a file
99 : : !!
100 : : !! OUTPUT
101 : : !! Only printing
102 : : !!
103 : : !! NOTES
104 : : !! bantot
105 : : !! doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy.
106 : : !! dom=frequency range
107 : : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
108 : : !! eigen11(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
109 : : !! in direction x
110 : : !! eigen12(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
111 : : !! in direction y
112 : : !! eigen13(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
113 : : !! in direction z
114 : : !! ecut=kinetic energy planewave cutoff (hartree).
115 : : !! fermie= fermi energy (Hartree)
116 : : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$).
117 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
118 : : !! kin11= Onsager kinetic coeficient=optical conductivity
119 : : !! kin12= Onsager kinetic coeficient
120 : : !! kin21= Onsager kinetic coeficient
121 : : !! kin22= Onsager kinetic coeficient
122 : : !! Kth=thermal conductivity
123 : : !! mom=number of frequency for conductivity computation
124 : : !! mband=maximum number of bands.
125 : : !! natom = number of atoms in the unit cell.
126 : : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
127 : : !! nkpt=number of k points in the IBZ for this perturbation
128 : : !! ngfft(3)=integer fft box dimensions.
129 : : !! nspinor=number of spinorial components of the wavefunctions.
130 : : !! nsppol=1 for unpolarized, 2 for spin-polarized.
131 : : !! ntypat = number of atom types.
132 : : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
133 : : !! occopt==option for occupancies
134 : : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor))
135 : : !! rprimd(3,3)=real space primitive translations.
136 : : !! of primitive translations.
137 : : !! Sth=thermopower
138 : : !! tsmear=smearing width (or temperature) in Hartree
139 : : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
140 : : !! wind=frequency windows for computations of sigma
141 : : !! wtk(nkpt)=weight assigned to each k point.
142 : : !! znucl(natom)=atomic number of atoms
143 : : !! np_sum=noziere-pines sumrule
144 : : !!
145 : : !! SOURCE
146 : :
147 : 14 : subroutine conducti_paw(filnam,filnam_out,varocc,use_gpu)
148 : :
149 : : !Arguments -----------------------------------
150 : : !scalars
151 : : character(len=fnlen),intent(in) :: filnam,filnam_out
152 : : integer, optional, intent(in) :: varocc
153 : : logical, optional, intent(in) :: use_gpu
154 : :
155 : : !Local variables-------------------------------
156 : : !scalars
157 : : integer,parameter :: master=0
158 : : integer :: bsize,bd_stride,dimid,iomode,bantot,bdtot_index,ncid,varid,nb_per_proc,etiq
159 : : integer :: comm,fform1,headform,iband,ijband,ierr,ikpt,master_band,idum,ib0,ib1,ij0,ij_stride
160 : : integer :: my_ncalls,ncalls_max,nb_blk
161 : : integer :: iom,isppol,jband,l1,l2,mband,me,mpierr,mom
162 : : integer :: natom,nband_k,nkpt,nproc,nspinor,nsppol,ntypat,broad_mode
163 : : integer :: occopt,iunt,opt_unt,occ_unt,iocc,pnp_size,add_drude,npairs,npairs_max,ipair,nkeep,nruns
164 : : integer :: lij_unt,sig_unt,sigd_unt,kth_unt,ocond_unt,occunit,occnpt,au_units
165 : : logical :: nc_unlimited,mykpt,myband,iomode_estf_mpiio,read_half_dipoles,omega0
166 : : logical :: my_use_gpu,run_sorted_k
167 : 14 : logical,allocatable :: band_keep(:),pair_drude(:),run_sorted(:)
168 : : real(dp) :: del,deltae,deltae_min,deltae_min_tmp,dhdk2_g,diff_eig,diff_occ
169 : : real(dp) :: dosdeltae,ecut,entropy,fermie,fermih,kin_fact,maxocc,docc_deig
170 : : real(dp) :: np_sum,np_sum_k1,np_sum_2,omin,omax,dom,oml,socc,socc_k,fact_omega0,pair_coef
171 : : real(dp) :: Tatm,tphysel,tsmear,ucvol,eig_in_max,eig_in_min,phi
172 : : character(len=fnlen) :: filnam1,filnam_gen,occfile
173 : : character(len=500) :: msg
174 : : character(len=100) :: line
175 : 14 : type(hdr_type) :: hdr
176 : 14 : type(MPI_type) :: mpi_enreg
177 : : !arrays
178 : : integer :: nc_count_5(5),nc_count_6(6),nc_start_5(5),nc_start_6(6),nc_stride_5(5),nc_stride_6(6)
179 : : integer :: dummy(0,0,0)
180 : 14 : integer,allocatable :: nband(:),run_start(:)
181 : : real(dp) :: dhdk2_r(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2)
182 : 14 : real(dp),allocatable :: cond_nd(:,:,:),cond_nd_k(:,:,:)
183 : 14 : real(dp),allocatable :: doccde(:),doccde_k(:),eig0_k(:),eigen0(:),eig0nc(:,:,:)
184 : 14 : real(dp),allocatable :: occ(:),occ_k(:),wtk(:),oml1(:),occ_in(:),eig_in(:),occ_tmp(:),ypp(:)
185 : 14 : real(dp),allocatable :: kin11(:,:),kin12(:),kin21(:),kin22(:)
186 : 14 : real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
187 : 14 : real(dp),allocatable :: psinablapsi(:,:,:),sig_abs(:),pair_ade(:),pair_wgt(:,:),band_np(:)
188 : :
189 : : ! *********************************************************************************
190 : :
191 : 14 : call timab(1186,1,tsec)
192 : : ABI_NVTX_START_RANGE(NVTX_CONDUCTI_PAW)
193 : :
194 [ + - ]: 14 : my_use_gpu=.false.; if (present(use_gpu)) my_use_gpu=use_gpu
195 : :
196 : : #if defined(HAVE_GPU_MARKERS)
197 : : call nvtxStartRange("read_input")
198 : : #endif
199 : : ! ---------------------------------------------------------------------------------
200 : : ! Read input data
201 : :
202 : : !Global MPI communicator
203 : 14 : comm = xmpi_world
204 : 14 : nproc = xmpi_comm_size(comm)
205 : 14 : me = xmpi_comm_rank(comm)
206 : :
207 : : !Read input parameters file
208 [ + + ]: 14 : if (me==master) then
209 [ - + ]: 8 : if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then
210 : 0 : ABI_ERROR(msg)
211 : : end if
212 : 8 : rewind(iunt)
213 : 8 : read(iunt,*)
214 : 8 : read(iunt,'(a)') filnam_gen ! Generic name for the files
215 [ + - ]: 8 : filnam1=trim(filnam_gen)//'_OPT'
216 : : ! Read frequency range
217 : 8 : read(iunt,*) dom,omin,omax,mom
218 [ + - ]: 8 : read(iunt,end=13,fmt='(a)') line
219 [ + + ]: 8 : read(line,end=11,fmt=*) broad_mode,au_units,phi,add_drude
220 : 1 : goto 14
221 [ - + ]: 7 : 11 read(line,end=13,fmt=*) broad_mode,au_units
222 : 0 : phi=zero; add_drude=0
223 : 0 : goto 14
224 : 7 : 13 broad_mode=1 ; au_units=0; phi=zero; add_drude=0
225 : : 14 continue
226 : : ! In case of varocc read filename of occupation datafile
227 [ - + ]: 8 : if (present(varocc)) then
228 [ # # ]: 0 : if (me==master) then
229 : 0 : write(std_out,*) 'Warning, this undocumented feature is highly experimental'
230 : 0 : write(std_out,*) 'and of limited physical validity, proceed with extreme caution!!!'
231 : : end if
232 : 0 : read(iunt,*) occfile !filename of the non-eq distribution function
233 : : end if
234 : 8 : close(iunt)
235 : : end if
236 : 14 : phi=phi*pi/180.0_dp
237 : :
238 : : !Send data to all procs
239 : 14 : call xmpi_bcast(dom,master,comm,mpierr)
240 : 14 : call xmpi_bcast(omin,master,comm,mpierr)
241 : 14 : call xmpi_bcast(omax,master,comm,mpierr)
242 : 14 : call xmpi_bcast(mom,master,comm,mpierr)
243 : 14 : call xmpi_bcast(broad_mode,master,comm,mpierr)
244 : 14 : call xmpi_bcast(au_units,master,comm,mpierr)
245 : 14 : call xmpi_bcast(phi,master,comm,mpierr)
246 : 14 : call xmpi_bcast(add_drude,master,comm,mpierr)
247 : :
248 : :
249 : : ! ---------------------------------------------------------------------------------
250 : : ! Read OPT file
251 : :
252 : : !Check for FORTRAN/.nc OPT file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF
253 [ + + ]: 14 : if (me==master) then
254 : 8 : call nctk_fort_or_ncfile(filnam1,iomode,msg)
255 [ - + ]: 8 : if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
256 : : end if
257 : 14 : call xmpi_bcast(filnam1,master,comm,mpierr)
258 : 14 : call xmpi_bcast(iomode,master,comm,mpierr)
259 : :
260 : : !Open OPT file and read HEADER
261 [ + + ]: 14 : if (me==master) then
262 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
263 [ - + ]: 8 : NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self))
264 : 8 : call hdr%ncread(ncid,fform1)
265 : : else
266 [ # # ]: 0 : if (open_file(filnam1,msg,newunit=opt_unt,form="unformatted",status="old")/=0) then
267 : 0 : ABI_ERROR(msg)
268 : : end if
269 : 0 : call hdr%fort_read(opt_unt,fform1,rewind=.true.)
270 : : end if
271 [ - + ]: 8 : ABI_CHECK(fform1/=0,sjoin("Error while reading ",filnam1))
272 [ - + ]: 8 : ABI_CHECK(fform1==610.or.fform1==620,"Conducti requires an OPT file with fform=610 or 620!")
273 : : end if
274 : 14 : call hdr%bcast(master,me,comm)
275 : 14 : call xmpi_bcast(fform1,master,comm,mpierr)
276 : :
277 : : !Extract info from the header
278 : 14 : headform=hdr%headform
279 : 14 : bantot=hdr%bantot
280 : 14 : ecut=hdr%ecut_eff
281 : 14 : natom=hdr%natom
282 : 14 : nkpt=hdr%nkpt
283 : 14 : nspinor=hdr%nspinor
284 : 14 : nsppol=hdr%nsppol
285 : 14 : ntypat=hdr%ntypat
286 : 14 : occopt=hdr%occopt
287 [ + + + + ]: 182 : rprimd(:,:)=hdr%rprimd(:,:)
288 : 14 : fermie=hdr%fermie
289 : 14 : tsmear=hdr%tsmear
290 [ - + + - : 42 : ABI_MALLOC(nband,(nkpt*nsppol))
+ - + - -
+ - + ]
291 [ - + + - : 42 : ABI_MALLOC(occ,(bantot))
+ - + - -
+ - + -
+ ]
292 [ - + + - : 42 : ABI_MALLOC(wtk,(nkpt))
+ - + - -
+ - + ]
293 [ + + ]: 30 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
294 [ + + ]: 248 : occ(1:bantot)=hdr%occ(1:bantot)
295 [ + + ]: 30 : wtk(1:nkpt)=hdr%wtk(1:nkpt)
296 [ + + ]: 30 : mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt)
297 : 14 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell
298 : :
299 : : !Read eigenvalues
300 [ - + + - : 42 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
+ - + - -
+ - + -
+ ]
301 [ + + ]: 14 : if (me==master) then
302 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
303 : 8 : varid=nctk_idname(ncid,"eigenvalues")
304 [ - + + - : 40 : ABI_MALLOC(eig0nc,(mband,nkpt,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
305 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,eig0nc))
306 [ + + + - : 16 : eigen0 = reshape(eig0nc,[mband*nkpt*nsppol])
+ - ]
307 [ - + ]: 8 : ABI_FREE(eig0nc)
308 : : !Close file here because the rest will possibly be read with collective I/O
309 [ - + ]: 8 : NCF_CHECK(nf90_close(ncid))
310 : : else
311 : 0 : read(opt_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
312 : : end if
313 : : end if
314 : 14 : call xmpi_bcast(eigen0,master,comm,mpierr)
315 : :
316 : : !In case of varocc, use arbitrary occupations
317 : : !Read occfile and overwrite occ by interpolating the data in the input file
318 [ - + ]: 14 : if (present(varocc)) then
319 [ # # ]: 0 : if (me==master) then
320 [ # # # # : 0 : ABI_MALLOC(occ_tmp,(bantot))
# # # # ]
321 [ # # ]: 0 : if (open_file(occfile,msg,newunit=occ_unt,form='formatted')/=0) then
322 : 0 : ABI_ERROR(msg)
323 : : end if
324 : : ! Get units used in occfile (1=Ha, 2= eV)
325 : 0 : read(occ_unt,*) occunit
326 : 0 : read(occ_unt,*) occnpt
327 [ # # # # : 0 : ABI_MALLOC(occ_in,(occnpt))
# # # # #
# # # ]
328 [ # # # # : 0 : ABI_MALLOC(eig_in,(occnpt))
# # # # ]
329 [ # # ]: 0 : do iocc=1,occnpt
330 : 0 : read(occ_unt,*) eig_in(iocc),occ_in(iocc)
331 [ # # ]: 0 : if (occunit==2) then !Convert units from eV to Hartree
332 : 0 : eig_in(iocc)=eig_in(iocc)/Ha_eV
333 : : occ_in(iocc)=occ_in(iocc)
334 : : end if
335 : : end do
336 : 0 : close(occ_unt)
337 : : ! Interpolation
338 [ # # # # : 0 : ABI_MALLOC(ypp,(occnpt))
# # # # #
# # # ]
339 : 0 : call spline(eig_in,occ_in,occnpt,zero,zero,ypp)
340 : : ! Interpolate neccessary values
341 [ # # # # : 0 : eig_in_max=maxval(eig_in)
# # # # ]
342 [ # # # # : 0 : eig_in_min=minval(eig_in)
# # # # ]
343 [ # # ]: 0 : do iocc=1,bantot
344 : : ! Check for Extrapolation and set them to physically sound values
345 [ # # ]: 0 : if (eigen0(iocc)<eig_in_min) then
346 : 0 : occ_tmp(iocc)=one
347 [ # # ]: 0 : else if (eigen0(iocc)>eig_in_max) then
348 : 0 : occ_tmp(iocc)=zero
349 : : else
350 : 0 : call splint(occnpt,eig_in,occ_in,ypp,1,eigen0(iocc),occ_tmp(iocc),ierr)
351 : : end if
352 : : end do
353 [ # # # # : 0 : occ=occ_tmp
# # # # ]
354 : : ! Clean up
355 [ # # ]: 0 : ABI_FREE(ypp)
356 [ # # ]: 0 : ABI_FREE(occ_tmp)
357 [ # # ]: 0 : ABI_FREE(occ_in)
358 [ # # ]: 0 : ABI_FREE(eig_in)
359 : : end if
360 : 0 : call xmpi_bcast(occ,master,comm,mpierr)
361 : 0 : occopt=2
362 : : end if ! varocc?
363 : :
364 : : #if defined(HAVE_GPU_MARKERS)
365 : : call nvtxEndRange()
366 : : #endif
367 : : #if defined(HAVE_GPU_MARKERS)
368 : : call nvtxStartRange("pre_spinkpt_loop")
369 : : #endif
370 : : !---------------------------------------------------------------------------------
371 : : ! Prepare kpt/band parallelization
372 : :
373 : 14 : call init_mpi_enreg(mpi_enreg)
374 : 14 : mpi_enreg%comm_kpt=comm
375 : 14 : mpi_enreg%me_kpt=me
376 : 14 : mpi_enreg%nproc_spkpt=nproc
377 : 14 : mpi_enreg%paralbd=1
378 [ - + + - : 70 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
379 [ - + + - : 42 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt))
+ - + - -
+ - + -
+ ]
380 : 14 : call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg)
381 : 14 : call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol)
382 : :
383 : : !---------------------------------------------------------------------------------
384 : : !Print some data
385 : :
386 : 14 : Tatm=tsmear*Ha_K
387 [ + + ]: 14 : if (me==master) then
388 : 8 : write(std_out,*)
389 : 8 : write(std_out,'(a)' )' Input data:'
390 : 8 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
391 : 8 : write(std_out,*)
392 : 8 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
393 : 8 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
394 : 8 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
395 : 8 : write(std_out,'(a,i8)') ' natom =',natom
396 : 8 : write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol
397 : 8 : write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
398 : 8 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
399 : 8 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
400 : : end if
401 : :
402 : : ! ---------------------------------------------------------------------------------
403 : : ! Compute derivative of occupations wrt the energy
404 : :
405 [ - + + - : 42 : ABI_MALLOC(doccde,(mband*nkpt*nsppol))
+ - + - -
+ - + ]
406 [ - + ]: 14 : if (occopt<=2) then
407 [ # # ]: 0 : if (me==master) then
408 : 0 : write(std_out,'(a,i4)') ' occopt =',occopt
409 : : end if
410 [ # # ]: 0 : doccde=zero
411 : : else
412 : 14 : tphysel=zero
413 : 14 : maxocc=two/(nsppol*nspinor)
414 : : call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermih,maxocc,mband,nband,&
415 : 14 : & socc,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,12,wtk)
416 : 14 : entropy=tsmear*entropy
417 [ + + ]: 14 : if (me==master) then
418 : 8 : write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha) =',entropy
419 : : end if
420 : 14 : entropy=entropy/socc
421 [ + + ]: 14 : if (me==master) then
422 : 8 : write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha/e) =',entropy
423 : : end if
424 : : endif
425 : :
426 : : !---------------------------------------------------------------------------------
427 : : ! Determine the frequency range and allocate frequency-dependent arrays
428 : :
429 : 14 : del=(omax-omin)/(mom-1)
430 [ - + + - : 42 : ABI_MALLOC(oml1,(mom))
+ - + - -
+ - + ]
431 [ + + ]: 14014 : do iom=1,mom
432 : 14014 : oml1(iom)=omin+dble(iom-1)*del
433 : : end do
434 : :
435 [ - + + - : 56 : ABI_MALLOC(kin11,(mom,nsppol))
- + + - +
- + - - +
- + - + ]
436 [ - + + - : 28 : ABI_MALLOC(kin12,(mom))
- + - + -
+ ]
437 [ - + + - : 28 : ABI_MALLOC(kin21,(mom))
- + - + -
+ ]
438 [ - + + - : 28 : ABI_MALLOC(kin22,(mom))
- + - + -
+ ]
439 [ - + + - : 42 : ABI_MALLOC(cond_nd,(3,3,mom))
+ - + - -
+ - + -
+ ]
440 [ - + + - : 28 : ABI_MALLOC(sig_abs,(mom))
- + - + ]
441 [ - + + - : 28 : ABI_MALLOC(Kth,(mom))
- + - + ]
442 [ - + + - : 28 : ABI_MALLOC(Stp,(mom))
- + - + ]
443 [ + + + + ]: 14028 : kin11 = zero
444 [ + + ]: 14014 : kin12 = zero
445 [ + + ]: 14014 : kin21 = zero
446 [ + + ]: 14014 : kin22 = zero
447 [ + + + + : 182014 : cond_nd = zero
+ + ]
448 [ + + ]: 14014 : sig_abs = zero
449 [ + + ]: 14014 : Kth = zero
450 [ + + ]: 14014 : Stp = zero
451 : :
452 : : #if defined(HAVE_GPU_MARKERS)
453 : : call nvtxStartRange("prepare_netcdf")
454 : : #endif
455 : : !---------------------------------------------------------------------------------
456 : : !Prepare valence-valence dipoles reading
457 : :
458 [ + - - + ]: 14 : iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
459 : :
460 : : !In case of netCDF access to OPT file, prepare collective I/O
461 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
462 [ + - ]: 14 : if (iomode_estf_mpiio) then
463 [ - + ]: 14 : NCF_CHECK(nctk_open_read(ncid,filnam1,comm))
464 : 14 : varid=nctk_idname(ncid,"dipole_valence_valence")
465 : 14 : nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR)
466 : 14 : read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR)
467 : : !Dipoles are read by blocks of (at most nb_blk) contiguous bands of the current proc
468 : : ! (nb_per_proc is the max number of bands of a k-point treated by a proc)
469 : 14 : nb_blk=int(dipoles_blk_mem/(48._dp*real(mband,dp)))
470 : 14 : nb_blk=max(1,min(nb_blk,nb_per_proc,mband))
471 [ + + ]: 14 : if (nproc>1) then
472 [ - + ]: 8 : NCF_CHECK(nctk_set_collective(ncid,varid))
473 : : !Count the reads of current proc, to complete them with dummy collective reads
474 : 8 : my_ncalls=count_band_blocks(mpi_enreg%proc_distrb,nband,nkpt,nsppol,me,nb_blk)
475 : 8 : call xmpi_max(my_ncalls,ncalls_max,comm,ierr)
476 : : end if
477 : : else
478 [ # # ]: 0 : if (me==master) then
479 [ # # ]: 0 : NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self))
480 : 0 : varid=nctk_idname(ncid,"dipole_valence_valence")
481 : : !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
482 : : ! NCF_CHECK(nctk_set_collective(ncid,varid))
483 : : !end if
484 : 0 : nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR)
485 : 0 : read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR)
486 : : end if
487 : 0 : call xmpi_bcast(nc_unlimited,master,comm,ierr)
488 : 0 : call xmpi_bcast(read_half_dipoles,master,comm,ierr)
489 : : end if
490 [ - + - - ]: 14 : if (nc_unlimited.and.read_half_dipoles) then
491 : 0 : msg="The OPT file has a wrong format!"
492 : 0 : ABI_BUG(msg)
493 : : end if
494 : : else
495 : 0 : read_half_dipoles=(fform1==620)
496 : : end if
497 : :
498 [ + - ]: 14 : if (iomode_estf_mpiio) then
499 : : !If MPI-IO, store only the elements of a block of nb_blk bands
500 [ - + + - : 42 : ABI_MALLOC(psinablapsi,(2,3,nb_blk*mband))
+ - + - -
+ - + -
+ ]
501 : : else
502 : : !If not, store all pairs (or half)
503 [ # # ]: 0 : if (read_half_dipoles) then ! only ib>=jb
504 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,(mband*(mband+1))/2))
# # # # #
# # # # #
# # ]
505 : : else
506 [ # # # # : 0 : ABI_MALLOC(psinablapsi,(2,3,mband*mband))
# # # # #
# # # #
# ]
507 : : end if
508 : : end if
509 [ + + ]: 56 : pnp_size=size(psinablapsi)
510 : :
511 : : #if defined(HAVE_GPU_MARKERS)
512 : : call nvtxEndRange()
513 : : #endif
514 : : !---------------------------------------------------------------------------------
515 : : ! Compute conductivity
516 : :
517 [ - + + - : 42 : ABI_MALLOC(kin11_k,(mom))
+ - + - -
+ - + ]
518 [ - + + - : 28 : ABI_MALLOC(kin12_k,(mom))
- + - + ]
519 [ - + + - : 28 : ABI_MALLOC(kin21_k,(mom))
- + - + ]
520 [ - + + - : 28 : ABI_MALLOC(kin22_k,(mom))
- + - + ]
521 [ - + + - : 42 : ABI_MALLOC(cond_nd_k,(3,3,mom))
+ - + - -
+ - + ]
522 : :
523 : : !Buffer of (n,m) band pairs: |e_n-e_m|, Drude flag and weights
524 : : ! (KG formula is applied to the buffer when full, see conducti_paw_kg_sum)
525 [ + - ]: 14 : npairs_max=max(merge(kg_pairs_blk_gpu,kg_pairs_blk,my_use_gpu),mband)
526 [ - + + - : 42 : ABI_MALLOC(pair_ade,(npairs_max))
+ - - + -
+ - + ]
527 [ - + + - : 42 : ABI_MALLOC(pair_drude,(npairs_max))
+ - - + -
+ - + ]
528 [ - + + - : 42 : ABI_MALLOC(pair_wgt,(9,npairs_max))
+ - - + -
+ - + ]
529 : : !Contributions to the sum rule and flags of the (n,m) pairs of band n
530 [ - + + - : 42 : ABI_MALLOC(band_np,(mband))
+ - + - -
+ - + ]
531 [ - + + - : 42 : ABI_MALLOC(band_keep,(mband))
+ - + - -
+ - + ]
532 : : !The buffer is made of runs: the pairs (n,m) of a band n, in the order of m.
533 : : ! As eigenvalues are sorted, |e_n-e_m| decreases along a run (checked in run_sorted)
534 [ + - - + : 42 : ABI_MALLOC(run_start,(mband+1))
+ - + - +
- - + -
+ ]
535 [ - + + - : 28 : ABI_MALLOC(run_sorted,(mband))
- + - + ]
536 : :
537 : 14 : np_sum = zero
538 : 14 : socc = zero
539 : 14 : deltae = zero
540 : 14 : deltae_min = 1.d99
541 : :
542 : : #if defined(HAVE_GPU_MARKERS)
543 : : call nvtxEndRange()
544 : : #endif
545 : : !LOOP OVER SPINS/K
546 : 14 : bdtot_index = 0
547 [ + + ]: 28 : do isppol=1,nsppol
548 [ + + ]: 44 : do ikpt=1,nkpt
549 : : #if defined(HAVE_GPU_MARKERS)
550 : : call nvtxStartRange("k-point step")
551 : : #endif
552 : 16 : etiq=ikpt+(isppol-1)*nkpt
553 : 16 : nband_k=nband(ikpt+(isppol-1)*nkpt)
554 : 16 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))
555 [ + + ]: 250 : master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
556 : :
557 : : #if defined(HAVE_GPU_MARKERS)
558 : : call nvtxStartRange("read_dipoles")
559 : : #endif
560 : : ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point
561 : : ! Master node reads and send to relevant processor
562 [ - + - - ]: 16 : if (.not.iomode_estf_mpiio.and.me==master) then
563 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
564 [ # # # # : 0 : psinablapsi=zero
# # ]
565 [ # # ]: 0 : if (nc_unlimited) then
566 [ # # # # ]: 0 : nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1]
567 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
568 [ # # ]: 0 : else if (.not.read_half_dipoles) then
569 [ # # # # ]: 0 : nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1]
570 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
571 : : else
572 [ # # # # ]: 0 : nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1]
573 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
574 : : end if
575 : : else
576 [ # # # # : 0 : psinablapsi=zero
# # ]
577 [ # # ]: 0 : bsize=nband_k**2;if (read_half_dipoles) bsize=(nband_k*(nband_k+1))/2
578 [ # # # # ]: 0 : read(opt_unt)(psinablapsi(1:2,1,ijband),ijband=1,bsize)
579 [ # # # # ]: 0 : read(opt_unt)(psinablapsi(1:2,2,ijband),ijband=1,bsize)
580 [ # # # # ]: 0 : read(opt_unt)(psinablapsi(1:2,3,ijband),ijband=1,bsize)
581 : : end if
582 [ # # ]: 0 : if (.not.mykpt) then
583 : 0 : call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,master_band,comm,etiq,ierr)
584 : : end if
585 : : end if
586 : : #if defined(HAVE_GPU_MARKERS)
587 : : call nvtxEndRange()
588 : : #endif
589 : :
590 : : ! Select k-points for current proc
591 [ + + ]: 16 : if (mykpt) then
592 : :
593 [ - + + - : 42 : ABI_MALLOC(eig0_k,(nband_k))
+ - + - -
+ - + ]
594 [ - + + - : 28 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
595 [ - + + - : 28 : ABI_MALLOC(doccde_k,(nband_k))
- + - + ]
596 : :
597 [ + + + + : 182014 : cond_nd_k = zero
+ + ]
598 [ + + ]: 14014 : kin11_k = zero
599 [ + + ]: 14014 : kin12_k = zero
600 [ + + ]: 14014 : kin21_k = zero
601 [ + + ]: 14014 : kin22_k = zero
602 : 14 : np_sum_k1 = zero
603 : 14 : socc_k = zero
604 : 14 : npairs = 0
605 : 14 : nruns = 0
606 : 14 : ib0 = 1
607 : 14 : ib1 = 0
608 : :
609 : : ! k-dependent data
610 [ + + ]: 228 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
611 [ + + ]: 228 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
612 [ + + ]: 228 : doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
613 : :
614 : : #if defined(HAVE_GPU_MARKERS)
615 : : call nvtxStartRange("dist_dipoles")
616 : : #endif
617 : : ! In case of non MPI-IO, receive all (n,m) dipoles from master proc
618 : : ! Then broadcast them to all band processors
619 [ - + ]: 14 : if (.not.iomode_estf_mpiio) then
620 [ # # # # ]: 0 : if (me/=master.and.me==master_band) then
621 : 0 : call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,me,comm,etiq,ierr)
622 : : end if
623 : 0 : call xmpi_bcast(psinablapsi,master,mpi_enreg%comm_band,mpierr)
624 : : end if
625 : : #if defined(HAVE_GPU_MARKERS)
626 : : call nvtxEndRange()
627 : : #endif
628 : :
629 : : #if defined(HAVE_GPU_MARKERS)
630 : : call nvtxStartRange("band_loop")
631 : : #endif
632 : : ! LOOP OVER BANDS n
633 [ + + ]: 228 : do iband=1,nband_k
634 : :
635 : : ! Select bands for current proc
636 : 214 : myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me)
637 [ + + ]: 228 : if (myband) then
638 : :
639 : : #if defined(HAVE_GPU_MARKERS)
640 : : call nvtxStartRange("read_dipoles_band")
641 : : #endif
642 : : ! In case of MPI-IO, read valence-valence dipoles for the block of bands ib0:ib1
643 : : ! starting at band n (contiguous bands of current proc)
644 [ + + ]: 144 : if (iomode_estf_mpiio.and.iband>ib1) then
645 : 14 : ib0=iband
646 : 14 : ib1=band_block_end(mpi_enreg%proc_distrb(ikpt,:,isppol),me,ib0,nband_k,nb_blk)
647 [ - + ]: 14 : if (nc_unlimited) then
648 [ # # # # ]: 0 : nc_start_6=[1,1,ib0,ikpt,isppol,1] ; nc_count_6=[2,3,ib1-ib0+1,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1]
649 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
650 [ - + ]: 14 : else if (.not.read_half_dipoles) then
651 [ # # # # ]: 0 : nc_start_6=[1,1,1,ib0,ikpt,isppol] ; nc_count_6=[2,3,mband,ib1-ib0+1,1,1] ; nc_stride_6=[1,1,1,1,1,1]
652 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
653 : : else
654 [ + + + + ]: 154 : nc_start_5=[1,1,(ib0*(ib0-1))/2+1,ikpt,isppol] ; nc_count_5=[2,3,(ib1*(ib1+1))/2-(ib0*(ib0-1))/2,1,1]
655 : 14 : nc_stride_5=[1,1,1,1,1]
656 [ - + ]: 14 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
657 : : end if
658 : : end if
659 : :
660 : : ! Position of the (n,m) dipoles in psinablapsi: ijband=ij0+(m-1)*ij_stride
661 : 144 : ij_stride=1
662 [ + - ]: 144 : if (read_half_dipoles) then
663 : : !Only n>=m pairs are stored; for MPI-IO, pairs of the block ib0:ib1
664 : 144 : ij0=(iband*(iband-1))/2-(ib0*(ib0-1))/2+1
665 [ # # # # ]: 0 : else if (iomode_estf_mpiio.and.nc_unlimited) then
666 : : !MPI-IO, unlimited bands dimension: psinablapsi(:,:,ib0:ib1,1:mband)
667 : 0 : ij0=iband-ib0+1
668 : 0 : ij_stride=ib1-ib0+1
669 : : else
670 : : !psinablapsi(:,:,1:bd_stride,ib0:ib1)
671 : : !bd_stride is mband for netCDF I/O, nband_k for Fortran I/O
672 [ # # ]: 0 : bd_stride=merge(mband,nband_k,iomode==IO_MODE_ETSF)
673 : 0 : ij0=(iband-ib0)*bd_stride+1
674 : : end if
675 : : #if defined(HAVE_GPU_MARKERS)
676 : : call nvtxEndRange()
677 : : #endif
678 : :
679 : : ! Apply KG formula to buffered pairs if there is no room for the pairs of band n
680 [ - + ]: 144 : if (npairs+iband>npairs_max) then
681 : 0 : run_start(nruns+1)=npairs+1
682 : : call conducti_paw_kg_sum(npairs,pair_ade,pair_drude,pair_wgt,nruns,run_start,run_sorted,&
683 : 0 : & mom,oml1,broad_mode,dom,my_use_gpu,kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
684 : 0 : npairs=0
685 : 0 : nruns=0
686 : : end if
687 : :
688 : : #if defined(HAVE_GPU_MARKERS)
689 : : call nvtxStartRange("fill_pairs")
690 : : #endif
691 : : ! LOOP OVER BANDS m
692 : : ! The (n,m) pair is stored in slot npairs+m of the buffer; pairs that do not
693 : : ! contribute are removed afterwards (serial loop, also computing deltae_min and
694 : : ! the sum rule in the order of m)
695 : : !$OMP PARALLEL DO IF(iband>=128) &
696 : : !$OMP PRIVATE(diff_occ,diff_eig,dhdk2_r,dhdk2_g,docc_deig,fact_omega0,ijband,ipair,kin_fact,l1,l2,omega0,pair_coef)
697 [ + + ]: 1328 : do jband=1,iband
698 : 1184 : diff_occ = occ_k(iband)-occ_k(jband)
699 : 1184 : diff_eig = eig0_k(iband)-eig0_k(jband)
700 : :
701 : 1184 : dhdk2_r = zero
702 : 1184 : dhdk2_g = zero
703 : :
704 : 1184 : ijband=ij0+(jband-1)*ij_stride
705 : :
706 [ + + ]: 4736 : do l2=1,3
707 [ + + ]: 15392 : do l1=1,3
708 : : dhdk2_r(l1,l2)=dhdk2_r(l1,l2)+(&
709 : : & psinablapsi(1,l1,ijband)*psinablapsi(1,l2,ijband)&
710 : 14208 : & +psinablapsi(2,l1,ijband)*psinablapsi(2,l2,ijband))
711 : : end do
712 : : end do
713 [ + + ]: 4736 : do l1=1,3
714 : : dhdk2_g=dhdk2_g &
715 : : & +(psinablapsi(1,l1,ijband)*psinablapsi(1,l1,ijband) &
716 : 4736 : & +psinablapsi(2,l1,ijband)*psinablapsi(2,l1,ijband))*third ! Average over directions
717 : : end do
718 : :
719 : : !Conductivity for each omega - Apply KG formula
720 : 1184 : kin_fact=(eig0_k(iband)+eig0_k(jband))*half-(fermie+entropy)
721 : 1184 : omega0=dabs(diff_occ)<tol12
722 [ + + ]: 1184 : fact_omega0=merge(two,one,omega0)
723 : : !Evaluate sumrule
724 : : if (.not.omega0) then
725 : 888 : docc_deig=dabs(diff_occ/diff_eig)
726 [ + + ]: 296 : else if(add_drude==1) then
727 : 13 : docc_deig=dabs(doccde_k(iband))
728 : : else
729 : : docc_deig=zero
730 : : endif
731 : 1184 : band_np(jband)=two*dhdk2_g*docc_deig/fact_omega0
732 : :
733 : : !Store the pair for KG formula
734 : : ! For omega/=0, docc_deig=|diff_occ|/omega: the 1/omega factor is applied in conducti_paw_kg_sum
735 : : ! Pairs with docc_deig=0 do not contribute
736 : 1184 : ipair=npairs+jband
737 : 1184 : pair_ade(ipair)=dabs(diff_eig)
738 [ + + + + ]: 1184 : band_keep(jband)=(.not.omega0).or.(add_drude==1)
739 [ + + ]: 1184 : if (.not.band_keep(jband)) cycle
740 [ + + ]: 901 : if (.not.omega0) then
741 : : pair_coef=dabs(diff_occ)
742 : : else
743 : 13 : pair_coef=dabs(doccde_k(iband))
744 : : endif
745 : 901 : pair_coef=pair_coef*pi/ucvol
746 : 901 : pair_drude(ipair)=omega0
747 : 901 : pair_wgt(1,ipair)=pair_coef*dhdk2_g
748 : 901 : pair_wgt(2,ipair)=pair_wgt(1,ipair)*kin_fact
749 : 901 : pair_wgt(3,ipair)=pair_wgt(2,ipair)*kin_fact
750 : 901 : pair_wgt(4,ipair)=pair_coef*dhdk2_r(1,1)
751 : 901 : pair_wgt(5,ipair)=pair_coef*dhdk2_r(2,2)
752 : 901 : pair_wgt(6,ipair)=pair_coef*dhdk2_r(3,3)
753 : 901 : pair_wgt(7,ipair)=pair_coef*dhdk2_r(1,2)
754 : 901 : pair_wgt(8,ipair)=pair_coef*dhdk2_r(1,3)
755 : 1328 : pair_wgt(9,ipair)=pair_coef*dhdk2_r(2,3)
756 : : end do !jband
757 : : !$OMP END PARALLEL DO
758 : :
759 : : ! Minimal validity limit, sum rule, and removal of the pairs that do not contribute
760 : : nkeep=0
761 : : run_sorted_k=.true.
762 [ + + ]: 1328 : do jband=1,iband
763 : 1184 : ipair=npairs+jband
764 : 1184 : deltae_min_tmp=pair_ade(ipair)
765 [ + + + + ]: 1184 : if ((deltae_min_tmp>=tol5).and.(deltae_min_tmp<=deltae_min)) deltae_min=deltae_min_tmp
766 : 1184 : np_sum_k1=np_sum_k1+band_np(jband)
767 [ + + ]: 1328 : if (band_keep(jband)) then
768 : 901 : nkeep=nkeep+1
769 [ + + ]: 901 : if (nkeep>1) then
770 [ + - ]: 786 : if (pair_ade(ipair)>pair_ade(npairs+nkeep-1)) run_sorted_k=.false.
771 : : end if
772 [ - + ]: 901 : if (nkeep<jband) then
773 : 0 : pair_ade(npairs+nkeep)=pair_ade(ipair)
774 : 0 : pair_drude(npairs+nkeep)=pair_drude(ipair)
775 [ # # ]: 0 : pair_wgt(:,npairs+nkeep)=pair_wgt(:,ipair)
776 : : end if
777 : : end if
778 : : end do
779 [ + + ]: 144 : if (nkeep>0) then
780 : 115 : nruns=nruns+1
781 : 115 : run_start(nruns)=npairs+1
782 : 115 : run_sorted(nruns)=run_sorted_k
783 : : end if
784 : 144 : npairs=npairs+nkeep
785 : : #if defined(HAVE_GPU_MARKERS)
786 : : call nvtxEndRange()
787 : : #endif
788 : :
789 : 144 : socc_k=socc_k+occ_k(iband)
790 : :
791 : : end if ! my band?
792 : : end do ! iband
793 : : #if defined(HAVE_GPU_MARKERS)
794 : : call nvtxEndRange()
795 : : #endif
796 : :
797 : : ! Apply KG formula to remaining buffered pairs
798 [ + - ]: 14 : if (npairs>0) then
799 : 14 : run_start(nruns+1)=npairs+1
800 : : call conducti_paw_kg_sum(npairs,pair_ade,pair_drude,pair_wgt,nruns,run_start,run_sorted,&
801 : 14 : & mom,oml1,broad_mode,dom,my_use_gpu,kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
802 : : end if
803 : :
804 : : ! Accumulate k-point contribution
805 [ + + ]: 14014 : do iom=1,mom
806 : 14000 : kin11(iom,isppol)=kin11(iom,isppol)+wtk(ikpt)*kin11_k(iom)
807 : 14000 : kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
808 : 14000 : kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
809 : 14000 : kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
810 [ + + + + ]: 182014 : cond_nd(:,:,iom)=cond_nd(:,:,iom)+wtk(ikpt)*cond_nd_k(:,:,iom)
811 : : end do
812 : 14 : np_sum=np_sum+wtk(ikpt)*np_sum_k1
813 : 14 : socc=socc+wtk(ikpt)*socc_k
814 : :
815 : : ! Validity limit
816 : 14 : deltae=deltae+(eig0_k(nband_k)-fermie)
817 : :
818 : 14 : ABI_FREE(eig0_k)
819 [ - + ]: 14 : ABI_FREE(occ_k)
820 [ - + ]: 14 : ABI_FREE(doccde_k)
821 : :
822 : : ! End loop over kpt/spin
823 : : end if ! My kpt?
824 : 30 : bdtot_index=bdtot_index+nband_k
825 : : #if defined(HAVE_GPU_MARKERS)
826 : : call nvtxEndRange()
827 : : #endif
828 : : end do ! ikpt
829 : : end do ! isppol
830 : :
831 : : #if defined(HAVE_GPU_MARKERS)
832 : : call nvtxStartRange("dummy_reads")
833 : : #endif
834 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
835 [ + - + + ]: 14 : if(iomode_estf_mpiio.and.nproc>1) then
836 [ + + ]: 10 : do idum=my_ncalls+1,ncalls_max
837 [ - + ]: 10 : if (nc_unlimited) then
838 : 0 : nc_start_6=[1,1,1,1,1,1] ; nc_count_6=[0,0,0,0,0,0] ; nc_stride_6=[1,1,1,1,1,1]
839 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
840 [ - + ]: 2 : else if (.not.read_half_dipoles) then
841 : 0 : nc_start_6=[1,1,1,1,1,1] ; nc_count_6=[0,0,0,0,0,0] ; nc_stride_6=[1,1,1,1,1,1]
842 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
843 : : else
844 : 2 : nc_start_5=[1,1,1,1,1] ; nc_count_5=[0,0,0,0,0] ; nc_stride_5=[1,1,1,1,1]
845 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
846 : : end if
847 : : enddo
848 : : endif
849 : : endif
850 : :
851 : : #if defined(HAVE_GPU_MARKERS)
852 : : call nvtxEndRange()
853 : : #endif
854 : :
855 : : #if defined(HAVE_GPU_MARKERS)
856 : : call nvtxStartRange("xmpi_sum")
857 : : #endif
858 : : !Accumulate kpt/band contributions over processors
859 : 14 : call xmpi_sum(kin11,comm,mpierr)
860 : 14 : call xmpi_sum(kin12,comm,mpierr)
861 : 14 : call xmpi_sum(kin21,comm,mpierr)
862 : 14 : call xmpi_sum(kin22,comm,mpierr)
863 : 14 : call xmpi_sum(np_sum,comm,mpierr)
864 : 14 : call xmpi_sum(socc,comm,mpierr)
865 : 14 : call xmpi_sum(deltae,comm,mpierr)
866 : 14 : call xmpi_min(deltae_min,comm,mpierr)
867 : 14 : call xmpi_sum(cond_nd,comm,mpierr)
868 : 14 : deltae=deltae/mpi_enreg%nproc_band
869 : : #if defined(HAVE_GPU_MARKERS)
870 : : call nvtxEndRange()
871 : : #endif
872 : :
873 [ - + ]: 14 : ABI_FREE(psinablapsi)
874 [ - + ]: 14 : ABI_FREE(pair_ade)
875 [ - + ]: 14 : ABI_FREE(pair_drude)
876 [ - + ]: 14 : ABI_FREE(pair_wgt)
877 [ - + ]: 14 : ABI_FREE(band_np)
878 [ - + ]: 14 : ABI_FREE(band_keep)
879 [ - + ]: 14 : ABI_FREE(run_start)
880 [ - + ]: 14 : ABI_FREE(run_sorted)
881 : :
882 : : #if defined(HAVE_GPU_MARKERS)
883 : : call nvtxStartRange("output")
884 : : #endif
885 : : !---------------------------------------------------------------------------------
886 : : ! Output results
887 : :
888 : : !Print file headers (only master node)
889 [ + + ]: 14 : if (me==master) then
890 : : ! Standard output
891 : : write(std_out,'(a,f10.5,a,f10.5,a)')&
892 : 8 : & ' Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha', &
893 : 16 : & deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
894 : : write(std_out,'(a,f10.5,a,f10.5,a)')&
895 : 8 : & ' DeltaE min =',deltae_min,' Ha',deltae_min*Ha_eV,' eV'
896 : :
897 : : ! _Lij file
898 [ - + ]: 8 : if (open_file(trim(filnam_out)//'_Lij',msg, newunit=lij_unt, form='formatted', action="write") /= 0) then
899 : 0 : ABI_ERROR(msg)
900 : : end if
901 : 8 : write(lij_unt,'(a)')' # omega(ua) L11 L12 L22'
902 : :
903 : : ! _sig_up_dn file
904 [ - + ]: 8 : if(nsppol==2) then
905 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_sig_up_dn', msg, newunit=sig_unt, form='formatted', action="write") /= 0) then
906 : 0 : ABI_ERROR(msg)
907 : : end if
908 [ # # ]: 0 : if(au_units>0) then
909 : 0 : write(sig_unt,'(2a)')' # omega(ua) cond(ua) ',&
910 : 0 : & ' cond(ua) UP cond(ua) DN'
911 : : else
912 : 0 : write(sig_unt,'(2a)')' # hbar*omega(eV) cond(ohm.cm)-1',&
913 : 0 : & ' cond(ohm.cm)-1 UP cond(ohm.cm)-1 DN'
914 : : endif
915 : : end if
916 : :
917 : : ! sig_dir file
918 [ - + ]: 8 : if (open_file(trim(filnam_out)//'_sig_tensor', msg, newunit=sigd_unt,form='formatted', action="write") /= 0) then
919 : 0 : ABI_ERROR(msg)
920 : : end if
921 [ - + ]: 8 : if(au_units>0) then
922 : 0 : write(sigd_unt,'(a)')' # omega(ua) cond_xx(ua) cond_yy(ua) cond_zz(ua) cond_xy(ua) cond_xz(ua) cond_yz(ua)'
923 : : else
924 : : write(sigd_unt,'(a)')' # hbar*omega(eV) cond_xx(ohm.cm)-1 cond_yy(ohm.cm)-1 cond_zz(ohm.cm)-1 &
925 : 8 : & cond_xy(ohm.cm)-1 cond_xz(ohm.cm)-1 cond_yz(ohm.cm)-1'
926 : : endif
927 : :
928 : : ! _Kth file
929 [ - + ]: 8 : if (open_file(trim(filnam_out)//'_Kth', msg, newunit=kth_unt, form='formatted', action="write") /=0) then
930 : 0 : ABI_ERROR(msg)
931 : : end if
932 : : write(kth_unt,'(a)')&
933 : 8 : & " #Thermal conductivity following B. Holst et al Phys. Rev. B 83 (2011) 235120"
934 [ - + ]: 8 : if(au_units>0) then
935 : 0 : write(kth_unt,'(a)')' # omega(ua) thermal cond(ua) thermopower(ua)'
936 : : else
937 : 8 : write(kth_unt,'(a)')' # hbar*omega(eV) thermal cond (W/m/K) thermopower(microohm/K)'
938 : : endif
939 : :
940 : : ! Output file
941 [ - + ]: 8 : if (open_file(trim(filnam_out)//'.out', msg, newunit=ocond_unt, form='formatted', action="write") /= 0) then
942 : 0 : ABI_ERROR(msg)
943 : : end if
944 : 8 : write(ocond_unt,'(a)' )'#Conducti output file:'
945 : 8 : write(ocond_unt,'(a)' )'#Contains all results produced by conducti utility'
946 : 8 : write(ocond_unt,'(a)' )'# '
947 : 8 : write(ocond_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =' ,mom,omin,omax,dom,' Ha'
948 : 8 : write(ocond_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1)
949 : 8 : write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2)
950 : 8 : write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3)
951 : 8 : write(ocond_unt,'(a,i8)' ) '# natom =',natom
952 : 8 : write(ocond_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol
953 : 8 : write(ocond_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha'
954 : 8 : write(ocond_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV'
955 : :
956 : 8 : write(ocond_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin'
957 : 8 : write(ocond_unt,'(a,f15.5)' )'# Number of electrons = ',socc
958 : 8 : write(ocond_unt,'(a,f15.5)' )'# sumrule = ',np_sum/socc
959 : : write(ocond_unt,'(a,f10.5,a,f10.5,a)' )&
960 : 8 : & '# Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
961 : 8 : write(ocond_unt,'(a)' )'# '
962 [ - + ]: 8 : if(au_units>0) then
963 : 0 : write(ocond_unt,'(a)')'# omega(ua) cond(ua) thermal cond(ua) thermopower(ua)'
964 : : else
965 : 8 : write(ocond_unt,'(a)')'# hbar*omega(eV) cond(ohm.cm)-1 thermal cond(W/m/K) thermopower(microvolt/K) '
966 : : endif
967 : :
968 : : end if ! me==master?
969 : :
970 : : !Compute (and print) thermal conductivity and thermopower
971 [ + + ]: 14014 : do iom=1,mom
972 : 14000 : oml=oml1(iom)
973 : :
974 [ + + ]: 28000 : do isppol=1,nsppol
975 [ + + ]: 14000 : if (dabs(kin11(iom,isppol))<tol19) kin11(iom,isppol)=zero
976 : 28000 : sig_abs(iom)=sig_abs(iom)+kin11(iom,isppol)
977 : : end do
978 : :
979 : 14000 : Kth(iom)=kin22(iom)/tsmear
980 : 14000 : Stp(iom)=zero
981 [ + + ]: 14000 : if (sig_abs(iom)/=zero) then
982 : 13169 : Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/sig_abs(iom))/tsmear
983 : 13169 : Stp(iom)=kin12(iom)/(sig_abs(iom)*tsmear)
984 : : else
985 : 831 : Kth(iom)=Kth(iom-1)
986 : 831 : Stp(iom)=Stp(iom-1)
987 : : endif
988 [ + + ]: 14000 : if (dabs(Kth(iom))<tol19) Kth(iom)=zero
989 [ - + ]: 14000 : if (dabs(Stp(iom))<tol19) Stp(iom)=zero
990 [ + + ]: 14000 : if (abs(kin12(iom))<10.0_dp**(-80)) kin12(iom)=zero
991 [ + + ]: 14000 : if (abs(kin21(iom))<10.0_dp**(-80)) kin21(iom)=zero
992 [ + + ]: 14014 : if (abs(kin22(iom))<10.0_dp**(-80)) kin22(iom)=zero
993 : : end do
994 : :
995 : : !Compute the imaginary part of the conductivity (principal value)
996 : : ! +derived optical properties.
997 [ + + ]: 14 : if (me==master) then
998 : 8 : call msig(sig_abs,mom,oml1,filnam_out,phi,au_units)
999 : : end if
1000 : :
1001 : : ! Units
1002 [ + + ]: 14 : if (me==master) then
1003 [ + + ]: 8008 : do iom=1,mom
1004 : 8000 : oml=oml1(iom)
1005 : 8008 : write(lij_unt,'(f12.5,3es22.12)') oml,sig_abs(iom),kin12(iom),kin22(iom)
1006 : : enddo
1007 : : endif
1008 : :
1009 : 14 : np_sum_2=zero
1010 [ + + ]: 14000 : do iom=1,mom-1
1011 : 14000 : np_sum_2=np_sum_2+(sig_abs(iom)+sig_abs(iom+1))*(oml1(iom+1)-oml1(iom))/2.0_dp
1012 : : enddo
1013 : 14 : np_sum_2=2.0_dp*np_sum_2*ucvol/pi
1014 [ + - + + ]: 14014 : if(au_units==0) oml1=oml1*Ha_eV
1015 [ + - + + : 182014 : if(au_units==0) cond_nd=cond_nd*Ohmcm
+ + + + ]
1016 [ + - + + ]: 14014 : if(au_units==0) sig_abs=sig_abs*Ohmcm
1017 [ + - + + : 14028 : if(au_units==0) kin11=kin11*Ohmcm
+ + ]
1018 [ + - + + ]: 14014 : if(au_units==0) Kth=Kth*3.4057d9/Ha_K
1019 [ + - + + ]: 14014 : if(au_units==0) Stp=Stp*Ha_J/(Ha_K*e_CB)*10.0_dp**6
1020 : :
1021 : :
1022 [ + + ]: 14 : if (me==master) then
1023 : 8 : write(std_out,'(a,f15.5)' )'# Number of electrons = ',socc
1024 : 8 : write(std_out,'(a,f15.5)' )'# sumrule = ',np_sum/socc
1025 : 8 : write(std_out,'(a,f15.5)' )'# sumrule (integration) = ',np_sum_2/socc
1026 [ + + ]: 8008 : do iom=1,mom
1027 : 8000 : oml=oml1(iom)
1028 : 8000 : write(sigd_unt,'(f12.5,6es22.12)') oml,cond_nd(1,1,iom),cond_nd(2,2,iom),cond_nd(3,3,iom),&
1029 : 16000 : & cond_nd(1,2,iom),cond_nd(1,3,iom),cond_nd(2,3,iom)
1030 [ - + ]: 8000 : if (nsppol==2) then
1031 : 0 : write(sig_unt,'(f12.5,3es22.12)') oml,sig_abs(iom), &
1032 : 0 : & kin11(iom,1),kin11(iom,2)
1033 : : end if
1034 : 8000 : write(kth_unt,'(f12.5,3es22.12)') oml,Kth(iom),Stp(iom)
1035 : 8008 : write(ocond_unt,'(1f12.5,3es22.12)') oml,sig_abs(iom),Kth(iom),Stp(iom)
1036 : : enddo
1037 : : end if
1038 : :
1039 : :
1040 : : !---------------------------------------------------------------------------------
1041 : : ! End
1042 : :
1043 : : !Close all files
1044 [ + + ]: 14 : if (me==master) then
1045 : 8 : write(std_out,'(2a)')ch10,'OUTPUT'
1046 : 8 : write(std_out,'(a)')trim(filnam_out)//'_Lij : Onsager kinetic coefficients'
1047 : 8 : write(std_out,'(a)')trim(filnam_out)//'_eps : Optical conductivity and dielectric function'
1048 : 8 : write(std_out,'(a)')trim(filnam_out)//'_sig_tensor : Optical conductivity tensor'
1049 : 8 : write(std_out,'(a)')trim(filnam_out)//'_Kth : Thermal conductivity and thermopower'
1050 : 8 : write(std_out,'(a)')trim(filnam_out)//'_abs : n, k, reflectivity, absorption'
1051 : 8 : close(lij_unt)
1052 [ - + ]: 8 : if(nsppol==2) close(sig_unt)
1053 : 8 : close(kth_unt)
1054 : 8 : close(ocond_unt)
1055 : 8 : close(sigd_unt)
1056 : : end if
1057 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
1058 [ - + - - ]: 14 : if (iomode_estf_mpiio.or.me==master) then
1059 [ - + ]: 14 : NCF_CHECK(nf90_close(ncid))
1060 : : end if
1061 [ # # ]: 0 : else if (me==master) then
1062 : 0 : ierr=close_unit(opt_unt,msg)
1063 [ # # ]: 0 : ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam1))
1064 : : end if
1065 : :
1066 : : #if defined(HAVE_GPU_MARKERS)
1067 : : call nvtxEndRange()
1068 : : #endif
1069 : : !Release memory space
1070 [ - + ]: 14 : ABI_FREE(kin11)
1071 [ - + ]: 14 : ABI_FREE(kin22)
1072 [ - + ]: 14 : ABI_FREE(kin12)
1073 [ - + ]: 14 : ABI_FREE(kin21)
1074 [ - + ]: 14 : ABI_FREE(kin11_k)
1075 [ - + ]: 14 : ABI_FREE(kin22_k)
1076 [ - + ]: 14 : ABI_FREE(kin12_k)
1077 [ - + ]: 14 : ABI_FREE(kin21_k)
1078 [ - + ]: 14 : ABI_FREE(Stp)
1079 [ - + ]: 14 : ABI_FREE(Kth)
1080 [ - + ]: 14 : ABI_FREE(cond_nd)
1081 [ - + ]: 14 : ABI_FREE(cond_nd_k)
1082 [ - + ]: 14 : ABI_FREE(sig_abs)
1083 [ - + ]: 14 : ABI_FREE(eigen0)
1084 [ - + ]: 14 : ABI_FREE(nband)
1085 [ - + ]: 14 : ABI_FREE(oml1)
1086 [ - + ]: 14 : ABI_FREE(occ)
1087 [ - + ]: 14 : ABI_FREE(doccde)
1088 [ - + ]: 14 : ABI_FREE(wtk)
1089 : 14 : call hdr%free()
1090 : 14 : call destroy_mpi_enreg(mpi_enreg)
1091 : :
1092 : : ABI_NVTX_END_RANGE()
1093 : 14 : call timab(1186,2,tsec)
1094 : :
1095 [ - + - + : 28 : end subroutine conducti_paw
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + ]
1096 : : !!***
1097 : :
1098 : : !----------------------------------------------------------------------
1099 : :
1100 : : !!****f* m_conducti/conducti_paw_kg_sum
1101 : : !! NAME
1102 : : !! conducti_paw_kg_sum
1103 : : !!
1104 : : !! FUNCTION
1105 : : !! Apply the Kubo-Greenwood formula to a buffer of (n,m) band pairs
1106 : : !! and accumulate the contributions for all frequencies.
1107 : : !! On CPU, parallelized over frequencies: each frequency sums over all pairs
1108 : : !! (gaussian broadening: only over the pairs of its window).
1109 : : !! On GPU (OpenMP offload), see conducti_paw_kg_sum_gpu.
1110 : : !!
1111 : : !! INPUTS
1112 : : !! npairs=number of pairs in the buffer
1113 : : !! pair_ade(npairs)=|e_n-e_m| for each pair
1114 : : !! pair_drude(npairs)=.true. for pairs with equal occupations (Drude term):
1115 : : !! their weights do not include the 1/omega factor
1116 : : !! pair_wgt(9,npairs)=weights of each pair (pi/ucvol*docc_deig already included,
1117 : : !! without the 1/omega factor):
1118 : : !! dhdk2_g, dhdk2_g*kin_fact, dhdk2_g*kin_fact**2,
1119 : : !! dhdk2_r(1,1),dhdk2_r(2,2),dhdk2_r(3,3),dhdk2_r(1,2),dhdk2_r(1,3),dhdk2_r(2,3)
1120 : : !! nruns=number of runs of pairs in the buffer
1121 : : !! run_start(nruns+1)=index of the first pair of each run (run_start(nruns+1)=npairs+1)
1122 : : !! run_sorted(nruns)=.true. if pair_ade decreases along the run
1123 : : !! mom=number of frequencies
1124 : : !! oml1(mom)=frequencies
1125 : : !! broad_mode=0 for gaussian broadening, 1 for lorentzian broadening
1126 : : !! dom=broadening width
1127 : : !! use_gpu=.true. to use GPU (OpenMP offload)
1128 : : !!
1129 : : !! SIDE EFFECTS
1130 : : !! kin11_k(mom),kin12_k(mom),kin21_k(mom),kin22_k(mom)=Onsager kinetic coefficients (accumulated)
1131 : : !! cond_nd_k(3,3,mom)=conductivity tensor (accumulated)
1132 : : !!
1133 : : !! SOURCE
1134 : :
1135 [ - + - + ]: 14 : subroutine conducti_paw_kg_sum(npairs,pair_ade,pair_drude,pair_wgt,nruns,run_start,run_sorted,&
1136 [ - + ]: 14 : & mom,oml1,broad_mode,dom,use_gpu,kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
1137 : :
1138 : : !Arguments -----------------------------------
1139 : : !scalars
1140 : : integer,intent(in) :: npairs,nruns,mom,broad_mode
1141 : : logical,intent(in) :: use_gpu
1142 : : real(dp),intent(in) :: dom
1143 : : !arrays
1144 : : integer,intent(in) :: run_start(nruns+1)
1145 : : logical,intent(in) :: pair_drude(:),run_sorted(nruns)
1146 : : real(dp),intent(in) :: pair_ade(:),pair_wgt(:,:),oml1(mom)
1147 : : real(dp),intent(inout) :: kin11_k(mom),kin12_k(mom),kin21_k(mom),kin22_k(mom),cond_nd_k(3,3,mom)
1148 : :
1149 : : !Local variables-------------------------------
1150 : : !scalars
1151 : : !exp(-x) is exactly zero in double precision for x>745.14: it is not evaluated beyond exp_arg_max
1152 : : ! (this avoids the slow underflow path of exp, without changing the results)
1153 : : real(dp),parameter :: exp_arg_max=746._dp
1154 : : integer :: iom,ipair,ipair1,ipair2,irun
1155 : : real(dp) :: arg1,arg2,arg_max,dirac,dirac_fact,dirac_width,inv_oml,oml,window
1156 : : !arrays
1157 : : real(dp) :: acc(9)
1158 : :
1159 : : ! *********************************************************************************
1160 : :
1161 [ + + ]: 14 : if (broad_mode==0) then
1162 : 1 : dirac_fact=one/(dom*sqrt(two_pi))
1163 : 1 : dirac_width=one/(sqrt(two)*dom)
1164 : : else
1165 : 13 : dirac_fact=dom/pi
1166 : 13 : dirac_width=dom**2
1167 : : end if
1168 : :
1169 : : #if defined(HAVE_GPU_MARKERS)
1170 : : call nvtxStartRange("kg_sum")
1171 : : #endif
1172 : :
1173 : : !Gaussian broadening: terms dirac_fact*exp(-arg) lower than 1d-20 are neglected (as in conducti_paw_core),
1174 : : ! i.e. terms with arg>arg_max (at most exp_arg_max, beyond which exp is exactly zero)
1175 : : ! For frequency omega, only pairs with ||e_n-e_m|-|omega||<=window contribute:
1176 : : ! beyond, both exp arguments are larger than arg_max
1177 : : ! In each run sorted by decreasing |e_n-e_m|, these pairs are found by binary search
1178 : 14 : arg_max=zero; window=zero
1179 [ + + ]: 14 : if (broad_mode==0) then
1180 : 1 : arg_max=min(exp_arg_max,max(zero,log(dirac_fact*1d20)))
1181 : 1 : window=sqrt(arg_max)/dirac_width
1182 : : end if
1183 : :
1184 : : #ifdef HAVE_OPENMP_OFFLOAD
1185 : : if (use_gpu) then
1186 : : call conducti_paw_kg_sum_gpu(npairs,pair_ade,pair_drude,pair_wgt,mom,oml1,broad_mode,&
1187 : : & dirac_fact,dirac_width,arg_max,kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
1188 : : else
1189 : : #else
1190 : : ABI_UNUSED(use_gpu)
1191 : : #endif
1192 : :
1193 : : !$OMP PARALLEL DO SCHEDULE(dynamic,8) PRIVATE(acc,arg1,arg2,dirac,inv_oml,ipair,ipair1,ipair2,irun,oml)
1194 [ + + ]: 14014 : do iom=1,mom
1195 : 14000 : oml=oml1(iom)
1196 : 14000 : inv_oml=one/abs(oml)
1197 : 14000 : acc(:)=zero
1198 [ + + ]: 14000 : if (broad_mode==0) then
1199 [ + + ]: 11000 : do irun=1,nruns
1200 : 10000 : ipair1=run_start(irun)
1201 : 10000 : ipair2=run_start(irun+1)-1
1202 [ + - ]: 10000 : if (run_sorted(irun)) then
1203 : 10000 : ipair1=run_start(irun)+count_sorted(pair_ade(run_start(irun):run_start(irun+1)-1),abs(oml)+window,.true.)
1204 : 10000 : ipair2=run_start(irun)-1+count_sorted(pair_ade(run_start(irun):run_start(irun+1)-1),abs(oml)-window,.false.)
1205 : : end if
1206 [ + + ]: 16358 : do ipair=ipair1,ipair2
1207 : : ! Take into account (n,m) and (m,n)
1208 : 5358 : arg1=((pair_ade(ipair)-oml)*dirac_width)**2
1209 : 5358 : arg2=((pair_ade(ipair)+oml)*dirac_width)**2
1210 [ - + - - ]: 5358 : if (arg1>arg_max.and.arg2>arg_max) cycle
1211 : 5358 : dirac=zero
1212 [ + - ]: 5358 : if (arg1<=arg_max) dirac=exp(-arg1)
1213 [ + + ]: 5358 : if (arg2<=arg_max) dirac=dirac-exp(-arg2)
1214 : 5358 : dirac=dirac_fact*dirac
1215 [ + + ]: 5358 : if (.not.pair_drude(ipair)) dirac=dirac*inv_oml
1216 [ + + ]: 63580 : acc(:)=acc(:)+dirac*pair_wgt(:,ipair)
1217 : : end do
1218 : : end do
1219 : : else
1220 [ + + ]: 859000 : do ipair=1,npairs
1221 : : dirac=dirac_fact*(one/((pair_ade(ipair)-oml)**2+dirac_width) &
1222 : 846000 : & -one/((pair_ade(ipair)+oml)**2+dirac_width))
1223 [ + - ]: 846000 : if (.not.pair_drude(ipair)) dirac=dirac*inv_oml
1224 [ + + ]: 8473000 : acc(:)=acc(:)+dirac*pair_wgt(:,ipair)
1225 : : end do
1226 : : end if
1227 : 14000 : kin11_k(iom)=kin11_k(iom)+acc(1)
1228 : 14000 : kin12_k(iom)=kin12_k(iom)-acc(2)
1229 : 14000 : kin21_k(iom)=kin21_k(iom)-acc(2)
1230 : 14000 : kin22_k(iom)=kin22_k(iom)+acc(3)
1231 : 14000 : cond_nd_k(1,1,iom)=cond_nd_k(1,1,iom)+acc(4)
1232 : 14000 : cond_nd_k(2,2,iom)=cond_nd_k(2,2,iom)+acc(5)
1233 : 14000 : cond_nd_k(3,3,iom)=cond_nd_k(3,3,iom)+acc(6)
1234 : 14000 : cond_nd_k(1,2,iom)=cond_nd_k(1,2,iom)+acc(7)
1235 : 14000 : cond_nd_k(2,1,iom)=cond_nd_k(2,1,iom)+acc(7)
1236 : 14000 : cond_nd_k(1,3,iom)=cond_nd_k(1,3,iom)+acc(8)
1237 : 14000 : cond_nd_k(3,1,iom)=cond_nd_k(3,1,iom)+acc(8)
1238 : 14000 : cond_nd_k(2,3,iom)=cond_nd_k(2,3,iom)+acc(9)
1239 : 14014 : cond_nd_k(3,2,iom)=cond_nd_k(3,2,iom)+acc(9)
1240 : : end do
1241 : : !$OMP END PARALLEL DO
1242 : : #ifdef HAVE_OPENMP_OFFLOAD
1243 : : end if
1244 : : #endif
1245 : :
1246 : : #if defined(HAVE_GPU_MARKERS)
1247 : : call nvtxEndRange()
1248 : : #endif
1249 : :
1250 : 14 : end subroutine conducti_paw_kg_sum
1251 : : !!***
1252 : :
1253 : : !----------------------------------------------------------------------
1254 : :
1255 : : !!****f* m_conducti/conducti_paw_kg_sum_gpu
1256 : : !! NAME
1257 : : !! conducti_paw_kg_sum_gpu
1258 : : !!
1259 : : !! FUNCTION
1260 : : !! GPU (OpenMP offload) version of conducti_paw_kg_sum.
1261 : : !! Parallelized over frequencies and chunks of pairs: each thread sums over a chunk
1262 : : !! of pairs (no window), the partial sums of the chunks are added afterwards
1263 : : !! (in the order of the chunks) in a second kernel.
1264 : : !! Same pairs and same neglected terms as on CPU (see conducti_paw_kg_sum).
1265 : : !!
1266 : : !! INPUTS
1267 : : !! npairs,pair_ade,pair_drude,pair_wgt,mom,oml1,broad_mode=see conducti_paw_kg_sum
1268 : : !! dirac_fact,dirac_width=factor and width of the broadening function
1269 : : !! arg_max=gaussian broadening: max argument of exp (terms beyond are neglected)
1270 : : !!
1271 : : !! SIDE EFFECTS
1272 : : !! kin11_k(mom),kin12_k(mom),kin21_k(mom),kin22_k(mom),cond_nd_k(3,3,mom)=see conducti_paw_kg_sum
1273 : : !!
1274 : : !! SOURCE
1275 : :
1276 : : #ifdef HAVE_OPENMP_OFFLOAD
1277 : : subroutine conducti_paw_kg_sum_gpu(npairs,pair_ade,pair_drude,pair_wgt,mom,oml1,broad_mode,&
1278 : : & dirac_fact,dirac_width,arg_max,kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
1279 : :
1280 : : !Arguments -----------------------------------
1281 : : !scalars
1282 : : integer,intent(in) :: npairs,mom,broad_mode
1283 : : real(dp),intent(in) :: dirac_fact,dirac_width,arg_max
1284 : : !arrays
1285 : : logical,intent(in) :: pair_drude(:)
1286 : : real(dp),intent(in) :: pair_ade(:),pair_wgt(:,:),oml1(mom)
1287 : : real(dp),intent(inout) :: kin11_k(mom),kin12_k(mom),kin21_k(mom),kin22_k(mom),cond_nd_k(3,3,mom)
1288 : :
1289 : : !Local variables-------------------------------
1290 : : !scalars
1291 : : integer :: chunk_size,ichunk,iom,ipair,nchunk
1292 : : real(dp) :: acc1,acc2,acc3,acc4,acc5,acc6,acc7,acc8,acc9,arg1,arg2,dirac,inv_oml,oml
1293 : : !arrays
1294 : : real(dp),allocatable :: partial(:,:,:)
1295 : :
1296 : : ! *********************************************************************************
1297 : :
1298 : : if (npairs==0) return
1299 : :
1300 : : !Chunks of pairs: partial(iom,:,ichunk) is the sum over the pairs of chunk ichunk
1301 : : chunk_size=(npairs+kg_gpu_nchunk-1)/kg_gpu_nchunk
1302 : : nchunk=(npairs+chunk_size-1)/chunk_size
1303 : : ABI_MALLOC(partial,(mom,9,nchunk))
1304 : : !$OMP TARGET DATA MAP(to:pair_ade(1:npairs),pair_drude(1:npairs),pair_wgt(:,1:npairs),oml1) &
1305 : : !$OMP& MAP(alloc:partial) MAP(tofrom:kin11_k,kin12_k,kin21_k,kin22_k,cond_nd_k)
1306 : : !Consecutive threads treat consecutive frequencies for the same chunk of pairs
1307 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) &
1308 : : !$OMP& PRIVATE(acc1,acc2,acc3,acc4,acc5,acc6,acc7,acc8,acc9,arg1,arg2,dirac,inv_oml,ipair,oml)
1309 : : do ichunk=1,nchunk
1310 : : do iom=1,mom
1311 : : oml=oml1(iom)
1312 : : inv_oml=one/abs(oml)
1313 : : acc1=zero; acc2=zero; acc3=zero; acc4=zero; acc5=zero; acc6=zero; acc7=zero; acc8=zero; acc9=zero
1314 : : do ipair=(ichunk-1)*chunk_size+1,min(ichunk*chunk_size,npairs)
1315 : : if (broad_mode==0) then
1316 : : ! Take into account (n,m) and (m,n)
1317 : : arg1=((pair_ade(ipair)-oml)*dirac_width)**2
1318 : : arg2=((pair_ade(ipair)+oml)*dirac_width)**2
1319 : : if (arg1>arg_max.and.arg2>arg_max) cycle
1320 : : dirac=zero
1321 : : if (arg1<=arg_max) dirac=exp(-arg1)
1322 : : if (arg2<=arg_max) dirac=dirac-exp(-arg2)
1323 : : dirac=dirac_fact*dirac
1324 : : else
1325 : : dirac=dirac_fact*(one/((pair_ade(ipair)-oml)**2+dirac_width) &
1326 : : & -one/((pair_ade(ipair)+oml)**2+dirac_width))
1327 : : end if
1328 : : if (.not.pair_drude(ipair)) dirac=dirac*inv_oml
1329 : : acc1=acc1+dirac*pair_wgt(1,ipair)
1330 : : acc2=acc2+dirac*pair_wgt(2,ipair)
1331 : : acc3=acc3+dirac*pair_wgt(3,ipair)
1332 : : acc4=acc4+dirac*pair_wgt(4,ipair)
1333 : : acc5=acc5+dirac*pair_wgt(5,ipair)
1334 : : acc6=acc6+dirac*pair_wgt(6,ipair)
1335 : : acc7=acc7+dirac*pair_wgt(7,ipair)
1336 : : acc8=acc8+dirac*pair_wgt(8,ipair)
1337 : : acc9=acc9+dirac*pair_wgt(9,ipair)
1338 : : end do
1339 : : partial(iom,1,ichunk)=acc1; partial(iom,2,ichunk)=acc2; partial(iom,3,ichunk)=acc3
1340 : : partial(iom,4,ichunk)=acc4; partial(iom,5,ichunk)=acc5; partial(iom,6,ichunk)=acc6
1341 : : partial(iom,7,ichunk)=acc7; partial(iom,8,ichunk)=acc8; partial(iom,9,ichunk)=acc9
1342 : : end do
1343 : : end do
1344 : : !Add the partial sums of the chunks (in the order of the chunks)
1345 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO &
1346 : : !$OMP& PRIVATE(acc1,acc2,acc3,acc4,acc5,acc6,acc7,acc8,acc9,ichunk)
1347 : : do iom=1,mom
1348 : : acc1=zero; acc2=zero; acc3=zero; acc4=zero; acc5=zero; acc6=zero; acc7=zero; acc8=zero; acc9=zero
1349 : : do ichunk=1,nchunk
1350 : : acc1=acc1+partial(iom,1,ichunk); acc2=acc2+partial(iom,2,ichunk); acc3=acc3+partial(iom,3,ichunk)
1351 : : acc4=acc4+partial(iom,4,ichunk); acc5=acc5+partial(iom,5,ichunk); acc6=acc6+partial(iom,6,ichunk)
1352 : : acc7=acc7+partial(iom,7,ichunk); acc8=acc8+partial(iom,8,ichunk); acc9=acc9+partial(iom,9,ichunk)
1353 : : end do
1354 : : kin11_k(iom)=kin11_k(iom)+acc1
1355 : : kin12_k(iom)=kin12_k(iom)-acc2
1356 : : kin21_k(iom)=kin21_k(iom)-acc2
1357 : : kin22_k(iom)=kin22_k(iom)+acc3
1358 : : cond_nd_k(1,1,iom)=cond_nd_k(1,1,iom)+acc4
1359 : : cond_nd_k(2,2,iom)=cond_nd_k(2,2,iom)+acc5
1360 : : cond_nd_k(3,3,iom)=cond_nd_k(3,3,iom)+acc6
1361 : : cond_nd_k(1,2,iom)=cond_nd_k(1,2,iom)+acc7
1362 : : cond_nd_k(2,1,iom)=cond_nd_k(2,1,iom)+acc7
1363 : : cond_nd_k(1,3,iom)=cond_nd_k(1,3,iom)+acc8
1364 : : cond_nd_k(3,1,iom)=cond_nd_k(3,1,iom)+acc8
1365 : : cond_nd_k(2,3,iom)=cond_nd_k(2,3,iom)+acc9
1366 : : cond_nd_k(3,2,iom)=cond_nd_k(3,2,iom)+acc9
1367 : : end do
1368 : : !$OMP END TARGET DATA
1369 : : ABI_FREE(partial)
1370 : :
1371 : : end subroutine conducti_paw_kg_sum_gpu
1372 : : #endif
1373 : : !!***
1374 : :
1375 : : !----------------------------------------------------------------------
1376 : :
1377 : : !!****f* m_conducti/band_block_end
1378 : : !! NAME
1379 : : !! band_block_end
1380 : : !!
1381 : : !! FUNCTION
1382 : : !! Return the last band of the block of bands starting at band ib0:
1383 : : !! contiguous bands treated by current proc, at most nb_blk bands.
1384 : : !!
1385 : : !! INPUTS
1386 : : !! distrb_k(:)=proc treating each band of the current k-point (proc_distrb)
1387 : : !! me=index of current proc
1388 : : !! ib0=first band of the block (treated by current proc)
1389 : : !! nband_k=number of bands of the current k-point
1390 : : !! nb_blk=maximum number of bands in a block
1391 : : !!
1392 : : !! SOURCE
1393 : :
1394 [ - + ]: 40 : pure integer function band_block_end(distrb_k,me,ib0,nband_k,nb_blk) result(ib1)
1395 : :
1396 : : !Arguments -----------------------------------
1397 : : integer,intent(in) :: me,ib0,nband_k,nb_blk
1398 : : integer,intent(in) :: distrb_k(:)
1399 : :
1400 : : ! *********************************************************************************
1401 : :
1402 : 40 : ib1=ib0
1403 [ + + + + ]: 358 : do while (ib1<nband_k.and.ib1-ib0+1<nb_blk)
1404 [ - + ]: 318 : if (distrb_k(ib1+1)/=me) exit
1405 : 40 : ib1=ib1+1
1406 : : end do
1407 : :
1408 : 40 : end function band_block_end
1409 : : !!***
1410 : :
1411 : : !----------------------------------------------------------------------
1412 : :
1413 : : !!****f* m_conducti/count_band_blocks
1414 : : !! NAME
1415 : : !! count_band_blocks
1416 : : !!
1417 : : !! FUNCTION
1418 : : !! Return the number of blocks of bands treated by current proc
1419 : : !! (see band_block_end), over all k-points and spins.
1420 : : !! Used to complete the collective reads of a proc with dummy reads.
1421 : : !!
1422 : : !! INPUTS
1423 : : !! proc_distrb(nkpt,mband,nsppol)=proc treating each band of each k-point and spin
1424 : : !! nband(nkpt*nsppol)=number of bands for each k-point and spin
1425 : : !! nkpt=number of k-points
1426 : : !! nsppol=number of spins
1427 : : !! me=index of current proc
1428 : : !! nb_blk=maximum number of bands in a block
1429 : : !!
1430 : : !! SOURCE
1431 : :
1432 [ - + ]: 16 : pure integer function count_band_blocks(proc_distrb,nband,nkpt,nsppol,me,nb_blk) result(nblocks)
1433 : :
1434 : : !Arguments -----------------------------------
1435 : : integer,intent(in) :: nkpt,nsppol,me,nb_blk
1436 : : integer,intent(in) :: proc_distrb(:,:,:),nband(nkpt*nsppol)
1437 : :
1438 : : !Local variables-------------------------------
1439 : : integer :: iband,ikpt,isppol,nband_k
1440 : :
1441 : : ! *********************************************************************************
1442 : :
1443 : 16 : nblocks=0
1444 [ + + ]: 32 : do isppol=1,nsppol
1445 [ + + ]: 48 : do ikpt=1,nkpt
1446 : 16 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1447 : 16 : iband=1
1448 [ + + ]: 224 : do while (iband<=nband_k)
1449 [ + + ]: 208 : if (proc_distrb(ikpt,iband,isppol)==me) then
1450 : 12 : nblocks=nblocks+1
1451 : 12 : iband=band_block_end(proc_distrb(ikpt,:,isppol),me,iband,nband_k,nb_blk)+1
1452 : : else
1453 : 180 : iband=iband+1
1454 : : end if
1455 : : end do
1456 : : end do
1457 : : end do
1458 : :
1459 : 16 : end function count_band_blocks
1460 : : !!***
1461 : :
1462 : : !----------------------------------------------------------------------
1463 : :
1464 : : !!****f* m_conducti/count_sorted
1465 : : !! NAME
1466 : : !! count_sorted
1467 : : !!
1468 : : !! FUNCTION
1469 : : !! Return the number of elements of an array sorted in descending order
1470 : : !! that are greater than x (strict=.true.) or greater or equal to x (strict=.false.).
1471 : : !!
1472 : : !! INPUTS
1473 : : !! list(:)=array sorted in descending order
1474 : : !! x=value
1475 : : !! strict=.true. to count elements >x, .false. to count elements >=x
1476 : : !!
1477 : : !! SOURCE
1478 : :
1479 [ - + ]: 20000 : pure integer function count_sorted(list,x,strict) result(nn)
1480 : :
1481 : : !Arguments -----------------------------------
1482 : : logical,intent(in) :: strict
1483 : : real(dp),intent(in) :: x
1484 : : real(dp),intent(in) :: list(:)
1485 : :
1486 : : !Local variables-------------------------------
1487 : : integer :: ihigh,ilow,imid
1488 : : logical :: above
1489 : :
1490 : : ! *********************************************************************************
1491 : :
1492 : : !Invariant: list(1:ilow) are counted, list(ihigh+1:) are not
1493 : 20000 : ilow=0
1494 : 20000 : ihigh=size(list)
1495 [ + + ]: 64825 : do while (ilow<ihigh)
1496 : 44825 : imid=(ilow+ihigh+1)/2
1497 [ + + ]: 44825 : if (strict) then
1498 : 22022 : above=(list(imid)>x)
1499 : : else
1500 : 22803 : above=(list(imid)>=x)
1501 : : end if
1502 [ + + ]: 64825 : if (above) then
1503 : : ilow=imid
1504 : : else
1505 : 41559 : ihigh=imid-1
1506 : : end if
1507 : : end do
1508 : 20000 : nn=ilow
1509 : :
1510 : 20000 : end function count_sorted
1511 : : !!***
1512 : :
1513 : : !----------------------------------------------------------------------
1514 : :
1515 : : !!****f* m_conducti/conducti_paw_core
1516 : : !! NAME
1517 : : !! conducti_paw_core
1518 : : !!
1519 : : !! FUNCTION
1520 : : !! This program computes the elements of the optical frequency dependent
1521 : : !! conductivity tensor and the conductivity along the three principal axes
1522 : : !! from the Kubo-Greenwood formula for PAW formalism
1523 : : !!
1524 : : !! INPUTS
1525 : : !! filnam=generic name for input data
1526 : : !! filnam_out=generic name for output data
1527 : : !! [with_absorption]=optiona lflag to activate the computation of absorption (_sigX file) (default=TRUE)
1528 : : !! [with_emissivity]=optiona lflag to activate the computation of emissivity (_emisX file) (default=FALSE)
1529 : : !!
1530 : : !! OUTPUT
1531 : : !! Only printing
1532 : : !!
1533 : : !! NOTES
1534 : : !! bantot
1535 : : !! dom=frequency range
1536 : : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
1537 : : !! ecut=kinetic energy planewave cutoff (hartree).
1538 : : !! fermie= fermi energy (Hartree)
1539 : : !! mom=number of frequency for conductivity computation
1540 : : !! mband=maximum number of bands.
1541 : : !! natom = number of atoms in the unit cell.
1542 : : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
1543 : : !! nkpt=number of k points in the IBZ for this perturbation
1544 : : !! ngfft(3)=integer fft box dimensions.
1545 : : !! nspinor=number of spinorial components of the wavefunctions.
1546 : : !! nsppol=1 for unpolarized, 2 for spin-polarized.
1547 : : !! ntypat = number of atom types.
1548 : : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
1549 : : !! occopt==option for occupancies
1550 : : !! psinablapsi2(2,3,mband,nphicor,natom)Matrix elements = <Phi_core|Nabla|Phi_i>
1551 : : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor))
1552 : : !! rprimd(3,3)=real space primitive translations.
1553 : : !! of primitive translations.
1554 : : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
1555 : : !! wind=frequency window for computations of sigma
1556 : : !! wtk(nkpt)=weight assigned to each k point.
1557 : : !!
1558 : : !! SOURCE
1559 : :
1560 : 14 : subroutine conducti_paw_core(filnam,filnam_out,with_absorption,with_emissivity)
1561 : :
1562 : : !Arguments -----------------------------------
1563 : : !scalars
1564 : : character(len=fnlen),intent(in) :: filnam,filnam_out
1565 : : logical,intent(in),optional :: with_absorption,with_emissivity
1566 : : !Local variables-------------------------------
1567 : : !scalars
1568 : : integer,parameter :: master=0
1569 : : integer :: iomode,atnbr,bantot,bdtot_index,comm,my_iband,input_atm,au_units
1570 : : integer :: ib0,ib1,my_ncalls,ncalls_max,nb_blk
1571 : : integer :: fform2,headform,iatom,iband,icor,ierr,ikpt,iatom_atnbr,itypat,itypat_atnbr
1572 : : integer :: iom,isppol,l1,mband,me,mom,mpierr,j2,etiq,pnp_size,broad_mode,absx_unt
1573 : : integer :: natom,nband_k,nkpt,nphicor,nproc,nspinor,nsppol,ntypat,nphicor_max,natom_atnbr
1574 : : integer :: occopt,iunt,opt2_unt,ncid,varid,master_band,nb_per_proc,idum,nbuf,nbuf_max
1575 : : integer :: sigx1_unt,sigx1_up_unt,sigx1_dn_unt,ems_unt,ems_up_unt,ems_dn_unt
1576 : : logical :: iomode_estf_mpiio,myband,mykpt,need_absorption,need_emissivity
1577 : : real(dp) :: del_sig,del_emis,deltae,ecut,fermie
1578 : : real(dp) :: omin,omax,omin_sig,omax_sig,omin_emis,omax_emis
1579 : : real(dp) :: dom,dom_ctr,dom_max,dom_tan1,dom_tan2,dhdk2_g
1580 : : real(dp) :: Tatm,tsmear,ucvol
1581 : : character(len=fnlen) :: filnam2,filnam_gen
1582 : : character(len=500) :: msg
1583 : : character(len=100) :: line
1584 : : character(len=20) :: str_atm
1585 : 14 : type(hdr_type) :: hdr
1586 : 14 : type(MPI_type) :: mpi_enreg
1587 : : !arrays
1588 : : integer :: nc_count(6),nc_start(6),nc_stride(6)
1589 : 14 : integer,allocatable :: nband(:),ncor(:,:),lcor(:,:),kappacor(:,:),typat(:),iatom_list(:)
1590 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),dummy(0,0,0,0,0),tsec(2)
1591 : 14 : real(dp),allocatable :: dom_var1(:,:),buf_eig(:),buf_wabs(:,:,:),buf_wemis(:,:,:),abs_diff_occ(:)
1592 : 14 : real(dp),allocatable :: eig0_k(:),eigen0(:),eig0nc(:,:,:)
1593 : 14 : real(dp),allocatable :: energy_cor(:,:),edge(:),nphicor_arr(:),occ_cor(:,:),maxocc_cor(:,:)
1594 : 14 : real(dp),allocatable :: occ(:),occ_k(:),wtk(:)
1595 : 14 : real(dp),allocatable :: oml_edge(:,:),oml_emis(:,:)
1596 : 14 : real(dp),allocatable :: psinablapsi2(:,:,:,:)
1597 : 14 : real(dp),allocatable :: sigx1(:,:,:,:),sigx1_av(:,:,:),sigx1_k(:,:,:)
1598 : 14 : real(dp),allocatable :: sum_spin_sigx1(:,:,:),sum_spin_sigx1_av(:,:)
1599 : 14 : real(dp),allocatable :: emisx(:,:,:,:),emisx_av(:,:,:),emisx_k(:,:,:)
1600 : 14 : real(dp),allocatable :: sum_spin_emisx(:,:,:),sum_spin_emisx_av(:,:)
1601 : :
1602 : : ! *********************************************************************************
1603 : : !optional flags
1604 [ + - ]: 14 : need_absorption=.true. ;if (present(with_absorption)) need_absorption=with_absorption
1605 [ + + ]: 14 : need_emissivity=.false.;if (present(with_emissivity)) need_emissivity=with_emissivity
1606 [ + - ]: 14 : if ((.not.need_absorption).and.(.not.need_emissivity)) return
1607 : :
1608 : 14 : call timab(1187,1,tsec)
1609 : : ABI_NVTX_START_RANGE(NVTX_CONDUCTI_PAW_CORE)
1610 : :
1611 : : #if defined(HAVE_GPU_MARKERS)
1612 : : call nvtxStartRange("core_read_input")
1613 : : #endif
1614 : : ! ---------------------------------------------------------------------------------
1615 : : ! Read input data
1616 : :
1617 : : !Global MPI communicators
1618 : 14 : comm = xmpi_world
1619 : 14 : nproc = xmpi_comm_size(comm)
1620 : 14 : me = xmpi_comm_rank(comm)
1621 : :
1622 : : !Read input parameters file
1623 [ + + ]: 14 : if (me==master) then
1624 [ - + ]: 8 : if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then
1625 : 0 : ABI_ERROR(msg)
1626 : : end if
1627 : 8 : rewind(iunt)
1628 : 8 : read(iunt,*)
1629 : 8 : read(iunt,'(a)')filnam_gen ! Generic name for the files
1630 [ + - ]: 8 : filnam2=trim(filnam_gen)//'_OPT2'
1631 : : ! Read frequency range
1632 [ + - ]: 8 : if (need_absorption) then
1633 : 8 : read(iunt,fmt='(a)') line
1634 [ + + ]: 8 : read(line,end=11,fmt=*) dom,omin,omax,mom,input_atm,dom_max,dom_ctr
1635 : 7 : goto 12
1636 : 1 : 11 read(line,fmt=*) dom,omin,omax,mom,input_atm ; dom_max=zero ; dom_ctr=zero
1637 : : 12 continue
1638 [ # # ]: 0 : else if (need_emissivity) then
1639 : 0 : read(iunt,*) dom,omin,omax,mom,input_atm
1640 : 0 : dom_max=zero;dom_ctr=zero
1641 : : end if
1642 [ + + ]: 8 : read(iunt,end=13,err=13,fmt=*) broad_mode,au_units
1643 : 1 : goto 14
1644 : 7 : 13 broad_mode=1 ; au_units=0
1645 : : 14 continue
1646 : 8 : close(iunt)
1647 [ + + - + ]: 8 : if (abs(dom_max)>tol10.and.dom_max<dom) then
1648 : 0 : msg = 'dom_max must be higher than dom!'
1649 : 0 : ABI_ERROR(msg)
1650 : : end if
1651 : : end if
1652 : :
1653 : : !Send data to all procs
1654 : 14 : call xmpi_bcast(dom,master,comm,mpierr)
1655 : 14 : call xmpi_bcast(omin,master,comm,mpierr)
1656 : 14 : call xmpi_bcast(omax,master,comm,mpierr)
1657 : 14 : call xmpi_bcast(mom,master,comm,mpierr)
1658 : 14 : call xmpi_bcast(input_atm,master,comm,mpierr)
1659 : 14 : call xmpi_bcast(dom_max,master,comm,mpierr)
1660 : 14 : call xmpi_bcast(dom_ctr,master,comm,mpierr)
1661 : 14 : call xmpi_bcast(broad_mode,master,comm,mpierr)
1662 : 14 : call xmpi_bcast(au_units,master,comm,mpierr)
1663 : :
1664 : :
1665 : : ! ---------------------------------------------------------------------------------
1666 : : ! Read OPT2 file
1667 : :
1668 : : !Check for FORTRAN/.nc OPT2 file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF
1669 [ + + ]: 14 : if (me==master) then
1670 : 8 : call nctk_fort_or_ncfile(filnam2,iomode,msg)
1671 [ - + ]: 8 : if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
1672 : : end if
1673 : 14 : call xmpi_bcast(filnam2,master,comm,mpierr)
1674 : 14 : call xmpi_bcast(iomode,master,comm,mpierr)
1675 : :
1676 : : !Open OPT2 file and read HEADER
1677 [ + + ]: 14 : if (me==master) then
1678 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
1679 [ - + ]: 8 : NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self))
1680 : 8 : call hdr%ncread(ncid,fform2)
1681 : : else
1682 [ # # ]: 0 : if (open_file(filnam2,msg,newunit=opt2_unt,form="unformatted",status="old")/=0) then
1683 : 0 : ABI_ERROR(msg)
1684 : : end if
1685 : 0 : call hdr%fort_read(opt2_unt,fform2,rewind=.true.)
1686 : : end if
1687 [ - + ]: 8 : ABI_CHECK(fform2/=0,sjoin("Error while reading ",filnam2))
1688 [ - + ]: 8 : ABI_CHECK(fform2==611.or.fform2==612.or.fform2==613,"OPT2 file format should be fform=611/612/613!")
1689 : : end if
1690 : 14 : call hdr%bcast(master,me,comm)
1691 : 14 : call xmpi_bcast(fform2,master,comm,mpierr)
1692 : :
1693 : : !Extract info from the header
1694 : 14 : headform=hdr%headform
1695 : 14 : bantot=hdr%bantot
1696 : 14 : ecut=hdr%ecut_eff
1697 : 14 : natom=hdr%natom
1698 : 14 : nkpt=hdr%nkpt
1699 : 14 : nspinor=hdr%nspinor
1700 : 14 : nsppol=hdr%nsppol
1701 : 14 : ntypat=hdr%ntypat
1702 [ - + + - : 42 : ABI_MALLOC(typat,(natom))
+ - + - -
+ - + ]
1703 [ + - - + : 77 : typat = hdr%typat
- - - - +
+ ]
1704 : : occopt=hdr%occopt
1705 [ + + + + ]: 182 : rprimd(:,:)=hdr%rprimd(:,:)
1706 : 14 : fermie=hdr%fermie
1707 : 14 : tsmear=hdr%tsmear
1708 [ - + + - : 42 : ABI_MALLOC(nband,(nkpt*nsppol))
+ - + - -
+ - + ]
1709 [ - + + - : 42 : ABI_MALLOC(occ,(bantot))
+ - + - -
+ - + ]
1710 [ - + + - : 42 : ABI_MALLOC(wtk,(nkpt))
+ - + - -
+ - + ]
1711 [ + + ]: 258 : occ(1:bantot)=hdr%occ(1:bantot)
1712 [ + + ]: 30 : wtk(1:nkpt)=hdr%wtk(1:nkpt)
1713 [ + + ]: 30 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1714 [ + + ]: 30 : mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt)
1715 : 14 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell
1716 : 14 : itypat_atnbr=typat(input_atm)
1717 : 14 : natom_atnbr=0
1718 [ + + ]: 63 : do iatom=1,natom
1719 [ + + ]: 49 : if(typat(iatom)==itypat_atnbr) natom_atnbr=natom_atnbr+1
1720 [ + + ]: 63 : if(iatom==input_atm) atnbr=natom_atnbr
1721 : : enddo
1722 : :
1723 : : !Read eigenvalues
1724 [ - + + - : 42 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
+ - + - -
+ - + -
+ ]
1725 [ + + ]: 14 : if (me==master) then
1726 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
1727 : 8 : varid=nctk_idname(ncid,"eigenvalues")
1728 [ - + + - : 40 : ABI_MALLOC(eig0nc,(mband,nkpt,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1729 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,eig0nc))
1730 [ + + + - : 16 : eigen0 = reshape(eig0nc,[mband*nkpt*nsppol])
+ - ]
1731 [ - + ]: 8 : ABI_FREE(eig0nc)
1732 : : else
1733 : 0 : read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
1734 : : end if
1735 : : end if
1736 : 14 : call xmpi_bcast(eigen0,master,comm,mpierr)
1737 : :
1738 : : !Read core states
1739 [ + + ]: 14 : if (me==master) then
1740 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
1741 [ - + ]: 8 : NCF_CHECK(nctk_get_dim(ncid,"max_number_of_core_states",nphicor_max))
1742 : : else
1743 : 0 : read(opt2_unt) nphicor_max
1744 : : end if
1745 : : end if
1746 : 14 : call xmpi_bcast(nphicor_max,master,comm,mpierr)
1747 : :
1748 [ - + + - : 56 : ABI_MALLOC(ncor,(nphicor_max,ntypat))
- + + - +
- + - - +
- + - + ]
1749 [ - + + - : 56 : ABI_MALLOC(occ_cor,(nphicor_max,ntypat))
- + + - +
- + - - +
- + - + ]
1750 [ - + + - : 42 : ABI_MALLOC(maxocc_cor,(nphicor_max,ntypat))
- + + - -
+ - + -
+ ]
1751 [ - + + - : 42 : ABI_MALLOC(lcor,(nphicor_max,ntypat))
- + + - -
+ - + -
+ ]
1752 [ - + + - : 42 : ABI_MALLOC(kappacor,(nphicor_max,ntypat))
- + + - -
+ - + -
+ ]
1753 [ - + + - : 42 : ABI_MALLOC(energy_cor,(nphicor_max,ntypat))
- + + - -
+ - + -
+ ]
1754 : :
1755 [ + + ]: 14 : if (me==master) then
1756 [ - + + - : 24 : ABI_MALLOC(nphicor_arr,(ntypat))
+ - + - -
+ - + -
+ ]
1757 [ + - ]: 8 : if (iomode==IO_MODE_ETSF) then
1758 : 8 : varid=nctk_idname(ncid,"n_quantum_number_core")
1759 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,ncor))
1760 : 8 : varid=nctk_idname(ncid,"l_quantum_number_core")
1761 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,lcor))
1762 : 8 : varid=nctk_idname(ncid,"kappa_core")
1763 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,kappacor))
1764 : 8 : varid=nctk_idname(ncid,"eigenvalues_core")
1765 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,energy_cor))
1766 : 8 : varid=nctk_idname(ncid,"occupation_core")
1767 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,occ_cor))
1768 : 8 : varid=nctk_idname(ncid,"max_occupation_core")
1769 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,maxocc_cor))
1770 : 8 : varid=nctk_idname(ncid,"number_of_core_states")
1771 [ - + ]: 8 : NCF_CHECK(nf90_get_var(ncid,varid,nphicor_arr))
1772 : : !Close here netcdf file here because the rest has to be read with collective I/O
1773 [ - + ]: 8 : NCF_CHECK(nf90_close(ncid))
1774 : : else
1775 [ # # ]: 0 : do itypat=1,ntypat
1776 : 0 : read(unit=opt2_unt) nphicor_arr(itypat)
1777 [ # # ]: 0 : do icor=1,nphicor_max
1778 : 0 : read(unit=opt2_unt,fmt='(a)') line
1779 : 0 : read(line,end=23,fmt=*) ncor(icor,itypat),lcor(icor,itypat),kappacor(icor,itypat),occ_cor(icor,itypat),maxocc_cor(icor,itypat),&
1780 [ # # ]: 0 : energy_cor(icor,itypat)
1781 : 0 : goto 24
1782 : 0 : 23 read(line,fmt=*) ncor(icor,itypat),lcor(icor,itypat),energy_cor(icor,itypat)
1783 : 0 : kappacor(icor,itypat)=0
1784 : 0 : occ_cor(icor,itypat)=one
1785 : 0 : maxocc_cor(icor,itypat)=one
1786 : 0 : 24 continue
1787 : : enddo
1788 : : end do
1789 : : end if
1790 : 8 : nphicor=nphicor_arr(itypat_atnbr)
1791 : 8 : ABI_FREE(nphicor_arr)
1792 : : end if ! master
1793 : 14 : call xmpi_bcast(ncor,master,comm,mpierr)
1794 : 14 : call xmpi_bcast(lcor,master,comm,mpierr)
1795 : 14 : call xmpi_bcast(kappacor,master,comm,mpierr)
1796 : 14 : call xmpi_bcast(energy_cor,master,comm,mpierr)
1797 : 14 : call xmpi_bcast(occ_cor,master,comm,mpierr)
1798 : 14 : call xmpi_bcast(maxocc_cor,master,comm,mpierr)
1799 : 14 : call xmpi_bcast(nphicor,master,comm,mpierr)
1800 : :
1801 [ - + + - : 42 : ABI_MALLOC(edge,(nphicor))
+ - + - -
+ - + ]
1802 [ + + ]: 62 : edge(1:nphicor)=fermie-energy_cor(1:nphicor,itypat_atnbr)
1803 : : #if defined(HAVE_GPU_MARKERS)
1804 : : call nvtxEndRange()
1805 : : #endif
1806 : : #if defined(HAVE_GPU_MARKERS)
1807 : : call nvtxStartRange("core_pre_spinkpt_loop")
1808 : : #endif
1809 : : !---------------------------------------------------------------------------------
1810 : : ! Prepare kpt/band parallelization
1811 : :
1812 : 14 : call init_mpi_enreg(mpi_enreg)
1813 : 14 : mpi_enreg%comm_kpt=comm
1814 : 14 : mpi_enreg%me_kpt=me
1815 : 14 : mpi_enreg%nproc_spkpt=nproc
1816 : 14 : mpi_enreg%paralbd=1
1817 [ - + + - : 70 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1818 [ - + + - : 42 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt))
+ - + - -
+ - + -
+ ]
1819 : 14 : call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg)
1820 : 14 : call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol)
1821 : : !---------------------------------------------------------------------------------
1822 : : !Print some data
1823 : 14 : Tatm=tsmear*Ha_K
1824 [ + + ]: 14 : if (me==master) then
1825 : 8 : write(std_out,*)
1826 : 8 : write(std_out,'(a)')'--------------------------------------------'
1827 : 8 : write(std_out,'(a,i4)') 'selected atom for X ray emission',input_atm
1828 : 8 : write(std_out,'(a)')'--------------------------------------------'
1829 [ + - ]: 8 : if (need_absorption) then
1830 [ + + ]: 8 : if (abs(dom_max)>tol10) then
1831 : 1 : write(std_out,'(a)')'************************ ARCTAN SMEARING'
1832 : 1 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1833 : 1 : write(std_out,'(a,2f10.5,a)')' dom_max,center =',dom_max,dom_ctr,' Ha'
1834 : : else
1835 : 7 : write(std_out,'(a)')'************************ FIXED SMEARING'
1836 : 7 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1837 : : endif
1838 : : end if
1839 : 8 : write(std_out,*)
1840 : 8 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1,1:3)
1841 : 8 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(2,1:3)
1842 : 8 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(3,1:3)
1843 : 8 : write(std_out,'(a,i8)') ' natom =',natom
1844 : 8 : write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol
1845 : 8 : write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
1846 : 8 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
1847 : 8 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
1848 : 8 : write(std_out,*)
1849 : 8 : write(std_out,*)
1850 : 8 : write(std_out,'(a)')'--------------------------------------------'
1851 : 8 : write(std_out,'(a,i4)') ' Number of core orbitals nc=',nphicor
1852 [ + + ]: 35 : do icor=1,nphicor
1853 [ + + ]: 35 : if (kappacor(icor,itypat_atnbr)==0) then
1854 : 15 : write(std_out,'(a,2i4,5f15.5)') ' n, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ', &
1855 : 15 : ncor(icor,itypat_atnbr),lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1856 : 15 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1857 : 30 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1858 : : else
1859 [ + + ]: 12 : if (kappacor(icor,itypat_atnbr)>0) then
1860 : 3 : j2=2*lcor(icor,itypat_atnbr)-1
1861 : 3 : write(std_out,'(a,i4,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1862 : 3 : & ncor(icor,itypat_atnbr),j2,' / 2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1863 : 3 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1864 : 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1865 : : else
1866 [ + + ]: 9 : if (kappacor(icor,itypat_atnbr)<-1) then
1867 : 3 : j2=2*lcor(icor,itypat_atnbr)+1
1868 : 3 : write(std_out,'(a,i4,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1869 : 3 : & ncor(icor,itypat_atnbr),j2,'/2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1870 : 3 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1871 : 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1872 : : else
1873 : 6 : write(std_out,'(a,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1874 : 6 : & ncor(icor,itypat_atnbr),' 1/2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1875 : 6 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1876 : 12 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1877 : : end if
1878 : : end if
1879 : : end if
1880 : : end do
1881 : 8 : write(std_out,'(a)')'--------------------------------------------'
1882 : : end if ! master
1883 : :
1884 : : !---------------------------------------------------------------------------------
1885 : : ! Determine the frequency range and allocate frequency-dependent arrays
1886 [ + - ]: 14 : if (need_absorption) then
1887 [ - + + - : 56 : ABI_MALLOC(oml_edge,(nphicor,mom))
- + + - +
- + - - +
- + ]
1888 [ - + + - : 42 : ABI_MALLOC(dom_var1,(nphicor,mom))
- + + - -
+ - + -
+ ]
1889 : 14 : omax_sig=omax ; omin_sig=omin
1890 : 14 : del_sig=(omax_sig-omin_sig)/(mom-1)
1891 [ + + ]: 14014 : do iom=1,mom
1892 [ + + ]: 62014 : do icor=1,nphicor
1893 : 48000 : oml_edge(icor,iom)=-energy_cor(icor,itypat_atnbr)+ dble(iom-1)*del_sig + omin_sig -one
1894 [ + + + + ]: 62000 : if ((oml_edge(icor,iom)<=edge(icor)).or.(abs(dom_max)<=tol10)) then
1895 : 45048 : dom_var1(icor,iom)= dom
1896 : : else
1897 : 2952 : dom_tan1= (oml_edge(icor,iom)-edge(icor))/dom_ctr
1898 : 2952 : dom_tan2=dom_tan1-one/dom_tan1**2
1899 : 2952 : dom_var1(icor,iom)=dom+dom_max*(half+piinv*datan(dom_tan2))
1900 : : endif
1901 : : enddo
1902 : : enddo
1903 [ - + + - : 84 : ABI_MALLOC(sigx1,(nphicor,mom,natom_atnbr,nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1904 [ + + + + : 209075 : sigx1=zero
+ + + + ]
1905 : : end if
1906 : :
1907 [ + + ]: 14 : if (need_emissivity) then
1908 [ - + + - : 48 : ABI_MALLOC(oml_emis,(nphicor,mom))
- + + - +
- + - - +
- + ]
1909 [ + - - + : 228 : omin_emis=minval(eigen0)
- - + + ]
1910 [ + - - + : 228 : omax_emis=maxval(eigen0)
- - + + ]
1911 : 12 : del_emis=(omax_emis-omin_emis)/(mom-1)
1912 [ + + ]: 12012 : do iom=1,mom
1913 [ + + ]: 52012 : do icor=1,nphicor
1914 : 52000 : oml_emis(icor,iom)=omin_emis+dble(iom-1)*del_emis!-energy_cor(icor,itypat_atnbr)
1915 : : enddo
1916 : : end do
1917 [ - + + - : 72 : ABI_MALLOC(emisx,(nphicor,mom,natom_atnbr,nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1918 [ + + + + : 184066 : emisx=zero
+ + + + ]
1919 : : end if
1920 : :
1921 : : #if defined(HAVE_GPU_MARKERS)
1922 : : call nvtxStartRange("core_prepare_netcdf")
1923 : : #endif
1924 : : !---------------------------------------------------------------------------------
1925 : : ! Prepare core-valence dipoles reading
1926 : :
1927 [ + - - + ]: 14 : iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
1928 : :
1929 : : !In case of netCDF access to OPT file, prepare collective I/O
1930 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
1931 [ + - ]: 14 : if (iomode_estf_mpiio) then
1932 [ - + ]: 14 : NCF_CHECK(nctk_open_read(ncid,filnam2,comm))
1933 : 14 : varid=nctk_idname(ncid,"dipole_core_valence")
1934 : : !Dipoles are read by blocks of (at most nb_blk) contiguous bands of the current proc
1935 : : ! (nb_per_proc is the max number of bands of a k-point treated by a proc)
1936 : 14 : nb_blk=int(dipoles_blk_mem/(24._dp*real(nphicor_max*natom,dp)))
1937 : 14 : nb_blk=max(1,min(nb_blk,nb_per_proc,mband))
1938 [ + + ]: 14 : if (nproc>1) then
1939 [ - + ]: 8 : NCF_CHECK(nctk_set_collective(ncid,varid))
1940 : : !Count the reads of current proc, to complete them with dummy collective reads
1941 : 8 : my_ncalls=count_band_blocks(mpi_enreg%proc_distrb,nband,nkpt,nsppol,me,nb_blk)
1942 : 8 : call xmpi_max(my_ncalls,ncalls_max,comm,ierr)
1943 : : end if
1944 [ # # ]: 0 : else if (me==master) then
1945 [ # # ]: 0 : NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self))
1946 : 0 : varid=nctk_idname(ncid,"dipole_core_valence")
1947 : : !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
1948 : : ! NCF_CHECK(nctk_set_collective(ncid,varid))
1949 : : !end if
1950 : : end if
1951 : : end if
1952 : :
1953 [ + - ]: 14 : if (iomode_estf_mpiio) then
1954 : : !If MPI-IO, store only the elements of a block of nb_blk bands
1955 [ - + + - : 70 : ABI_MALLOC(psinablapsi2,(3,nphicor_max,natom,nb_blk))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1956 : : else
1957 : : !If not, store the elements for all bands
1958 [ # # # # : 0 : ABI_MALLOC(psinablapsi2,(3,nphicor_max,natom,mband))
# # # # #
# # # # #
# # # # #
# # # ]
1959 : : end if
1960 [ + + ]: 70 : pnp_size=size(psinablapsi2)
1961 : : #if defined(HAVE_GPU_MARKERS)
1962 : : call nvtxEndRange()
1963 : : #endif
1964 : : !---------------------------------------------------------------------------------
1965 : : ! Compute X absorption coefficient and/or X emissivity
1966 : :
1967 [ + - ]: 14 : if (need_absorption) then
1968 [ - + + - : 70 : ABI_MALLOC(sigx1_k,(nphicor,mom,natom_atnbr))
- + + - -
+ + - + -
+ - - + -
+ ]
1969 : : end if
1970 [ + + ]: 14 : if (need_emissivity) then
1971 [ - + + - : 60 : ABI_MALLOC(emisx_k,(nphicor,mom,natom_atnbr))
- + + - -
+ + - + -
+ - - + -
+ ]
1972 : : end if
1973 : :
1974 : : !List of the atoms of the selected type
1975 [ - + + - : 42 : ABI_MALLOC(iatom_list,(natom_atnbr))
+ - + - -
+ - + ]
1976 : 14 : iatom_atnbr=0
1977 [ + + ]: 63 : do iatom=1,natom
1978 [ + + ]: 63 : if (typat(iatom)==itypat_atnbr) then
1979 : 47 : iatom_atnbr=iatom_atnbr+1
1980 : 47 : iatom_list(iatom_atnbr)=iatom
1981 : : end if
1982 : : end do
1983 : :
1984 : : !Buffer of bands: eigenvalue and weights for each core state and atom
1985 : : ! (weights include pi/ucvol and occupations, not the 1/omega factor)
1986 : : ! Contributions are computed when the buffer is full, see conducti_paw_core_sum
1987 : 14 : nbuf_max=min(core_bands_blk,mband)
1988 [ - + + - : 42 : ABI_MALLOC(buf_eig,(nbuf_max))
+ - + - -
+ - + -
+ ]
1989 [ + - ]: 14 : if (need_absorption) then
1990 [ - + + - : 70 : ABI_MALLOC(buf_wabs,(nbuf_max,nphicor,natom_atnbr))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1991 : : end if
1992 [ + + ]: 14 : if (need_emissivity) then
1993 [ - + + - : 60 : ABI_MALLOC(buf_wemis,(nbuf_max,nphicor,natom_atnbr))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1994 : : end if
1995 [ - + + - : 42 : ABI_MALLOC(abs_diff_occ,(nphicor))
+ - + - -
+ - + ]
1996 : :
1997 : 14 : deltae = zero
1998 : :
1999 : : #if defined(HAVE_GPU_MARKERS)
2000 : : call nvtxEndRange()
2001 : : #endif
2002 : : !LOOP OVER SPINS/K
2003 : 14 : bdtot_index = 0
2004 [ + + ]: 28 : do isppol=1,nsppol
2005 [ + + ]: 44 : do ikpt=1,nkpt
2006 : : #if defined(HAVE_GPU_MARKERS)
2007 : : call nvtxStartRange("core_k-point step")
2008 : : #endif
2009 : 16 : etiq=ikpt+(isppol-1)*nkpt
2010 : 16 : nband_k=nband(ikpt+(isppol-1)*nkpt)
2011 : 16 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))
2012 [ + + ]: 260 : master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
2013 : : #if defined(HAVE_GPU_MARKERS)
2014 : : call nvtxStartRange("core_read_dipoles")
2015 : : #endif
2016 : : ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point
2017 : : ! Master node reads and send to relevant processor
2018 [ - + - - ]: 16 : if (.not.iomode_estf_mpiio.and.me==master) then
2019 [ # # ]: 0 : if (iomode==IO_MODE_ETSF) then
2020 [ # # ]: 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
2021 [ # # ]: 0 : nc_count=[3,nphicor_max,natom,mband,1,1]
2022 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
2023 : : else
2024 [ # # # # : 0 : psinablapsi2=zero
# # # # ]
2025 [ # # ]: 0 : if (fform2==612) then ! New OPT2 file format
2026 : 0 : read(opt2_unt) (((psinablapsi2(1,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
2027 : 0 : read(opt2_unt) (((psinablapsi2(2,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
2028 : 0 : read(opt2_unt) (((psinablapsi2(3,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
2029 [ # # ]: 0 : else if (fform2==613) then ! Large OPT2 file format
2030 [ # # ]: 0 : do iband=1,nband_k
2031 : 0 : read(opt2_unt) ((psinablapsi2(1,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
2032 : 0 : read(opt2_unt) ((psinablapsi2(2,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
2033 : 0 : read(opt2_unt) ((psinablapsi2(3,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
2034 : : end do
2035 : : else
2036 : : !The old writing was not efficient (indexes order is bad)
2037 [ # # ]: 0 : do iatom=1,natom
2038 [ # # # # ]: 0 : read(opt2_unt) ((psinablapsi2(1,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
2039 [ # # # # ]: 0 : read(opt2_unt) ((psinablapsi2(2,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
2040 [ # # # # ]: 0 : read(opt2_unt) ((psinablapsi2(3,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
2041 : : end do
2042 : : end if
2043 : : end if
2044 [ # # ]: 0 : if (.not.mykpt) then
2045 : 0 : call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,master_band,comm,etiq,ierr)
2046 : : end if
2047 : : end if
2048 : : #if defined(HAVE_GPU_MARKERS)
2049 : : call nvtxEndRange()
2050 : : #endif
2051 : : !! Select k-points for current proc
2052 [ + + ]: 16 : if (mykpt) then
2053 : :
2054 [ - + + - : 42 : ABI_MALLOC(eig0_k,(nband_k))
+ - + - -
+ - + ]
2055 [ - + + - : 28 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
2056 : :
2057 [ + - + + : 185055 : if (need_absorption) sigx1_k=zero
+ + + + ]
2058 [ + + + + : 160050 : if (need_emissivity) emisx_k=zero
+ + + + ]
2059 : 14 : nbuf=0
2060 : 14 : ib0=1
2061 : 14 : ib1=0
2062 : :
2063 : : ! k-dependent data
2064 [ + + ]: 238 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
2065 [ + + ]: 238 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
2066 : :
2067 : : #if defined(HAVE_GPU_MARKERS)
2068 : : call nvtxStartRange("core_dist_dipoles")
2069 : : #endif
2070 : : ! In case of non MPI-IO, receive all (n,m) dipoles from master proc
2071 : : ! Then broadcast them to all band processors
2072 [ - + ]: 14 : if (.not.iomode_estf_mpiio) then
2073 [ # # # # ]: 0 : if (me/=master.and.me==master_band) then
2074 : 0 : call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,me,comm,etiq,ierr)
2075 : : end if
2076 : 0 : call xmpi_bcast(psinablapsi2,master,mpi_enreg%comm_band,mpierr)
2077 : : end if
2078 : : #if defined(HAVE_GPU_MARKERS)
2079 : : call nvtxEndRange()
2080 : : #endif
2081 : :
2082 : : #if defined(HAVE_GPU_MARKERS)
2083 : : call nvtxStartRange("core_band_loop")
2084 : : #endif
2085 : : ! LOOP OVER BANDS
2086 : :
2087 [ + + ]: 238 : do iband=1,nband_k
2088 : :
2089 : : ! Select bands for current proc
2090 : 224 : myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me)
2091 [ + + ]: 238 : if (myband) then
2092 : :
2093 : : #if defined(HAVE_GPU_MARKERS)
2094 : : call nvtxStartRange("core_read_dipoles_band")
2095 : : #endif
2096 : : ! In case of MPI-IO, read core-valence dipoles for the block of bands ib0:ib1
2097 : : ! starting at band n (contiguous bands of current proc)
2098 [ + + ]: 154 : if (iomode_estf_mpiio.and.iband>ib1) then
2099 : 14 : ib0=iband
2100 : 14 : ib1=band_block_end(mpi_enreg%proc_distrb(ikpt,:,isppol),me,ib0,nband_k,nb_blk)
2101 [ + + ]: 98 : nc_start=[1,1,1,ib0,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
2102 [ + + ]: 98 : nc_count=[3,nphicor_max,natom,ib1-ib0+1,1,1]
2103 [ - + ]: 14 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
2104 : : end if
2105 : : ! Position of band n in psinablapsi2: block ib0:ib1 for MPI-IO, all bands otherwise
2106 : 154 : my_iband=iband-ib0+1
2107 : : #if defined(HAVE_GPU_MARKERS)
2108 : : call nvtxEndRange()
2109 : : #endif
2110 : :
2111 : : ! Compute contributions of buffered bands if the buffer is full
2112 [ - + ]: 154 : if (nbuf==nbuf_max) then
2113 [ # # ]: 0 : if (need_absorption) then
2114 : : call conducti_paw_core_sum(nbuf,buf_eig,buf_wabs,nphicor,natom_atnbr,mom,energy_cor(:,itypat_atnbr),&
2115 : 0 : & oml_edge,broad_mode,dom,sigx1_k,lor_width=dom_var1)
2116 : : end if
2117 [ # # ]: 0 : if (need_emissivity) then
2118 : : call conducti_paw_core_sum(nbuf,buf_eig,buf_wemis,nphicor,natom_atnbr,mom,energy_cor(:,itypat_atnbr),&
2119 : 0 : & oml_edge,broad_mode,dom,emisx_k)
2120 : : end if
2121 : 0 : nbuf=0
2122 : : end if
2123 : :
2124 : : ! Store band n in the buffer
2125 : : ! For absorption, docc_deig=|diff_occ|/omega, for emission docc_deig=|occ_n|/omega:
2126 : : ! the 1/omega factor is applied in conducti_paw_core_sum
2127 : 154 : nbuf=nbuf+1
2128 : 154 : buf_eig(nbuf)=eig0_k(iband)
2129 [ + + ]: 676 : do icor=1,nphicor
2130 : : abs_diff_occ(icor)=abs((occ_cor(icor,itypat_atnbr)/maxocc_cor(icor,itypat_atnbr) &
2131 : 676 : & -occ_k(iband)/two*nsppol*nspinor)*two/nspinor/nsppol)
2132 : : end do
2133 [ + + ]: 518 : do iatom_atnbr=1,natom_atnbr
2134 : 364 : iatom=iatom_list(iatom_atnbr)
2135 [ + + ]: 1790 : do icor=1,nphicor
2136 : : dhdk2_g=zero
2137 [ + + ]: 5088 : do l1=1,3
2138 : 5088 : dhdk2_g=dhdk2_g+psinablapsi2(l1,icor,iatom,my_iband)*third
2139 : : end do
2140 [ + - ]: 1272 : if (need_absorption) buf_wabs(nbuf,icor,iatom_atnbr)=dhdk2_g*abs_diff_occ(icor)*pi/ucvol
2141 [ + + ]: 1636 : if (need_emissivity) buf_wemis(nbuf,icor,iatom_atnbr)=dhdk2_g*abs(occ_k(iband))*pi/ucvol
2142 : : end do
2143 : : end do
2144 : : end if ! my band?
2145 : : end do ! iband
2146 : : #if defined(HAVE_GPU_MARKERS)
2147 : : call nvtxEndRange()
2148 : : #endif
2149 : :
2150 : : ! Compute contributions of remaining buffered bands
2151 [ + - ]: 14 : if (nbuf>0) then
2152 [ + - ]: 14 : if (need_absorption) then
2153 : : call conducti_paw_core_sum(nbuf,buf_eig,buf_wabs,nphicor,natom_atnbr,mom,energy_cor(:,itypat_atnbr),&
2154 : 14 : & oml_edge,broad_mode,dom,sigx1_k,lor_width=dom_var1)
2155 : : end if
2156 [ + + ]: 14 : if (need_emissivity) then
2157 : : call conducti_paw_core_sum(nbuf,buf_eig,buf_wemis,nphicor,natom_atnbr,mom,energy_cor(:,itypat_atnbr),&
2158 : 12 : & oml_edge,broad_mode,dom,emisx_k)
2159 : : end if
2160 : : end if
2161 : :
2162 : : ! Accumulate k-point contribution
2163 [ + - ]: 14 : if (need_absorption) then
2164 : : sigx1(1:nphicor,1:mom,1:natom_atnbr,isppol)=sigx1(1:nphicor,1:mom,1:natom_atnbr,isppol) &
2165 [ + + + + : 185055 : & +wtk(ikpt)*sigx1_k(1:nphicor,1:mom,1:natom_atnbr)
+ + ]
2166 : : end if
2167 [ + + ]: 14 : if (need_emissivity) then
2168 : : emisx(1:nphicor,1:mom,1:natom_atnbr,isppol)=emisx(1:nphicor,1:mom,1:natom_atnbr,isppol) &
2169 [ + + + + : 160048 : & +wtk(ikpt)*emisx_k(1:nphicor,1:mom,1:natom_atnbr)
+ + ]
2170 : : end if
2171 : :
2172 : : ! Validity limit
2173 : 14 : deltae=deltae+eig0_k(nband_k)
2174 : :
2175 : 14 : ABI_FREE(eig0_k)
2176 [ - + ]: 14 : ABI_FREE(occ_k)
2177 : : ! End loop over kpt/spin
2178 : : end if ! My kpt?
2179 : 30 : bdtot_index=bdtot_index+nband_k
2180 : : #if defined(HAVE_GPU_MARKERS)
2181 : : call nvtxEndRange()
2182 : : #endif
2183 : : end do ! ikpt
2184 : : end do ! isppol
2185 : :
2186 : : #if defined(HAVE_GPU_MARKERS)
2187 : : call nvtxStartRange("core_dummy_reads")
2188 : : #endif
2189 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
2190 [ + - + + ]: 14 : if(iomode_estf_mpiio.and.nproc>1) then
2191 [ + + ]: 10 : do idum=my_ncalls+1,ncalls_max
2192 : 2 : nc_start=[1,1,1,1,1,1];nc_stride=[1,1,1,1,1,1]
2193 : 2 : nc_count=[0,0,0,0,0,0]
2194 [ - + ]: 10 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start,stride=nc_stride,count=nc_count))
2195 : : enddo
2196 : : endif
2197 : : endif
2198 : :
2199 : : #if defined(HAVE_GPU_MARKERS)
2200 : : call nvtxEndRange()
2201 : : #endif
2202 : :
2203 : : #if defined(HAVE_GPU_MARKERS)
2204 : : call nvtxStartRange("core_xmpi_sum")
2205 : : #endif
2206 : : !Accumulate kpt/band contributions over processors
2207 [ + - ]: 14 : if (need_absorption) then
2208 : 14 : call xmpi_sum(sigx1,comm,mpierr)
2209 : : end if
2210 [ + + ]: 14 : if (need_emissivity) then
2211 : 12 : call xmpi_sum(emisx,comm,mpierr)
2212 : : end if
2213 : 14 : call xmpi_sum(deltae,comm,mpierr)
2214 : 14 : deltae=deltae/mpi_enreg%nproc_band
2215 : : #if defined(HAVE_GPU_MARKERS)
2216 : : call nvtxEndRange()
2217 : : #endif
2218 : :
2219 : :
2220 : : #if defined(HAVE_GPU_MARKERS)
2221 : : call nvtxStartRange("core_post_process")
2222 : : #endif
2223 : : !Release some memory
2224 [ - + ]: 14 : ABI_FREE(iatom_list)
2225 [ - + ]: 14 : ABI_FREE(buf_eig)
2226 [ + - ]: 14 : if (need_absorption) then
2227 [ - + ]: 14 : ABI_FREE(buf_wabs)
2228 : : end if
2229 [ + + ]: 14 : if (need_emissivity) then
2230 [ - + ]: 12 : ABI_FREE(buf_wemis)
2231 : : end if
2232 [ - + ]: 14 : ABI_FREE(abs_diff_occ)
2233 [ - + ]: 14 : ABI_FREE(psinablapsi2)
2234 [ + - ]: 14 : if (need_absorption) then
2235 [ - + ]: 14 : ABI_FREE(sigx1_k)
2236 : : end if
2237 [ + + ]: 14 : if (need_emissivity) then
2238 [ - + ]: 12 : ABI_FREE(emisx_k)
2239 : : end if
2240 : :
2241 : : !Close core-valence dipoles file
2242 [ + - ]: 14 : if (iomode == IO_MODE_ETSF) then
2243 [ - + - - ]: 14 : if (iomode_estf_mpiio.or.me==master) then
2244 [ - + ]: 14 : NCF_CHECK(nf90_close(ncid))
2245 : : end if
2246 [ # # ]: 0 : else if (me==master) then
2247 : 0 : ierr=close_unit(opt2_unt,msg)
2248 [ # # ]: 0 : ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam2))
2249 : : end if
2250 : :
2251 : : !---------------------------------------------------------------------------------
2252 : : ! Post-processing
2253 : :
2254 [ + - ]: 14 : if (need_absorption) then
2255 [ + + + + : 209075 : sigx1=sigx1*dble(natom_atnbr)
+ + + + ]
2256 [ - + + - : 70 : ABI_MALLOC(sigx1_av,(nphicor,mom,nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
2257 [ + + + + : 62028 : sigx1_av=zero
+ + ]
2258 : : endif
2259 [ + + ]: 14 : if (need_emissivity) then
2260 : : ! Filter low values
2261 [ + + ]: 24 : do isppol=1,nsppol
2262 [ + + ]: 66 : do iatom_atnbr=1,natom_atnbr
2263 [ + + ]: 42054 : do iom=1,mom
2264 [ + + ]: 184042 : do icor=1,nphicor
2265 [ - + ]: 184000 : if (emisx(icor,iom,iatom_atnbr,isppol)<=tol16) emisx(icor,iom,iatom_atnbr,isppol)=zero
2266 : : end do
2267 : : end do
2268 : : end do
2269 : : end do
2270 [ + + + + : 184066 : emisx=emisx*dble(natom_atnbr)
+ + + + ]
2271 [ - + + - : 60 : ABI_MALLOC(emisx_av,(nphicor,mom,nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
2272 [ + + + + : 52024 : emisx_av=zero
+ + ]
2273 : : endif
2274 : :
2275 [ + + ]: 28 : do isppol=1,nsppol
2276 [ + + ]: 75 : do iatom_atnbr=1,natom_atnbr
2277 [ + + ]: 47061 : do iom=1,mom
2278 [ + + ]: 209047 : do icor=1,nphicor
2279 [ + - ]: 162000 : if(need_absorption) sigx1_av(icor,iom,isppol)=sigx1_av(icor,iom,isppol)+sigx1(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
2280 [ + + ]: 209000 : if(need_emissivity) emisx_av(icor,iom,isppol)=emisx_av(icor,iom,isppol)+emisx(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
2281 : : end do
2282 : : end do
2283 : : end do
2284 : : end do
2285 : :
2286 : : ! Spin treatment
2287 [ - + ]: 14 : if(nsppol==2) then
2288 [ # # ]: 0 : if(need_absorption) then
2289 [ # # # # : 0 : ABI_MALLOC(sum_spin_sigx1,(nphicor,mom,natom_atnbr))
# # # # #
# # # # #
# # # # #
# ]
2290 [ # # # # : 0 : ABI_MALLOC(sum_spin_sigx1_av,(nphicor,mom))
# # # # #
# # # # #
# # ]
2291 [ # # # # : 0 : sum_spin_sigx1=zero ; sum_spin_sigx1_av=zero
# # # # #
# ]
2292 : : endif
2293 [ # # ]: 0 : if(need_emissivity) then
2294 [ # # # # : 0 : ABI_MALLOC(sum_spin_emisx,(nphicor,mom,natom_atnbr))
# # # # #
# # # # #
# # # # #
# ]
2295 [ # # # # : 0 : ABI_MALLOC(sum_spin_emisx_av,(nphicor,mom))
# # # # #
# # # # #
# # ]
2296 [ # # # # : 0 : sum_spin_emisx=zero ; sum_spin_emisx_av=zero
# # # # #
# ]
2297 : : endif
2298 [ # # ]: 0 : do isppol=1,nsppol
2299 [ # # ]: 0 : do iatom_atnbr=1,natom_atnbr
2300 [ # # ]: 0 : do iom=1,mom
2301 [ # # ]: 0 : do icor=1,nphicor
2302 [ # # ]: 0 : if(need_absorption) sum_spin_sigx1(icor,iom,iatom_atnbr)=sum_spin_sigx1(icor,iom,iatom_atnbr) &
2303 : 0 : & +sigx1(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
2304 [ # # ]: 0 : if(need_emissivity) sum_spin_emisx(icor,iom,iatom_atnbr)=sum_spin_emisx(icor,iom,iatom_atnbr) &
2305 : 0 : & +emisx(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
2306 : : end do
2307 : : end do
2308 : : end do
2309 : : end do
2310 [ # # ]: 0 : do isppol=1,nsppol
2311 [ # # ]: 0 : do iom=1,mom
2312 [ # # ]: 0 : do icor=1,nphicor
2313 [ # # ]: 0 : if(need_absorption) sum_spin_sigx1_av(icor,iom)=sum_spin_sigx1_av(icor,iom)+sigx1_av(icor,iom,isppol)
2314 [ # # ]: 0 : if(need_emissivity) sum_spin_emisx_av(icor,iom)=sum_spin_emisx_av(icor,iom)+emisx_av(icor,iom,isppol)
2315 : : end do
2316 : : end do
2317 : : end do
2318 : : endif
2319 : :
2320 : : ! Units
2321 [ + - ]: 14 : if(au_units==0) then
2322 [ + + + + ]: 62014 : oml_edge=oml_edge*Ha_eV
2323 [ + - ]: 14 : if(need_absorption) then
2324 [ + + + + : 62028 : sigx1_av=sigx1_av*Ohmcm
+ + ]
2325 [ + + + + : 209075 : sigx1=sigx1*Ohmcm
+ + + + ]
2326 [ - + ]: 14 : if(nsppol==2) then
2327 [ # # # # ]: 0 : sum_spin_sigx1_av=sum_spin_sigx1_av*Ohmcm
2328 [ # # # # : 0 : sum_spin_sigx1=sum_spin_sigx1*Ohmcm
# # ]
2329 : : endif
2330 : : endif
2331 [ + + ]: 14 : if(need_emissivity) then
2332 [ + + + + : 52024 : emisx_av=emisx_av*Ohmcm
+ + ]
2333 [ + + + + : 184066 : emisx=emisx*Ohmcm
+ + + + ]
2334 [ - + ]: 12 : if(nsppol==2) then
2335 [ # # # # ]: 0 : sum_spin_emisx_av=sum_spin_emisx_av*Ohmcm
2336 [ # # # # : 0 : sum_spin_emisx=sum_spin_emisx*Ohmcm
# # ]
2337 : : endif
2338 : : endif
2339 : : endif
2340 : :
2341 : : #if defined(HAVE_GPU_MARKERS)
2342 : : call nvtxEndRange()
2343 : : #endif
2344 : : #if defined(HAVE_GPU_MARKERS)
2345 : : call nvtxStartRange("core_output")
2346 : : #endif
2347 : : !---------------------------------------------------------------------------------
2348 : : ! Output results
2349 : 14 : write(str_atm,*) input_atm
2350 : 14 : str_atm=adjustl(str_atm)
2351 : : str_atm=trim(str_atm)
2352 : :
2353 : : !Only master node outputs results in files (only master node)
2354 [ + + ]: 14 : if (me==master) then
2355 : :
2356 : : ! Standard output
2357 [ + - ]: 8 : if (need_absorption) then
2358 : 8 : write(std_out,*) 'Absorption: valence state orbital energies: omin,omax',omin_sig,omax_sig
2359 : : end if
2360 [ + + ]: 8 : if (need_emissivity) then
2361 : 6 : write(std_out,*) 'Emissivity: valence state orbital energies: omin,omax',omin_emis,omax_emis
2362 : : end if
2363 [ + - ]: 8 : if (need_absorption) then
2364 : : write(std_out,'(a,f10.5,a,f10.5,a)')&
2365 : 8 : & ' Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
2366 : : end if
2367 : :
2368 : : ! _sigX file
2369 : : if (need_absorption) then
2370 [ - + ]: 8 : if (open_file(trim(filnam_out)//'_sigX_at'//str_atm,msg,newunit=sigx1_unt,form='formatted',action="write")/=0) then
2371 : 0 : ABI_ERROR(msg)
2372 : : end if
2373 [ + + ]: 8 : if (abs(dom_max)>tol10) then
2374 : 1 : write(sigx1_unt,'(a)')'#***************************************************** ARCTAN SMEARING ********'
2375 : 1 : write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
2376 : 1 : write(sigx1_unt,'(a,2f10.5,a)')'# dom_max,center =',dom_max,dom_ctr,' Ha'
2377 : : else
2378 : 7 : write(sigx1_unt,'(a)')'#***************************************************** FIXED SMEARING ********'
2379 : 7 : write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
2380 : : endif
2381 : 8 : write(sigx1_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1)
2382 : 8 : write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2)
2383 : 8 : write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3)
2384 : 8 : write(sigx1_unt,'(a,i8)' ) '# natom =',natom
2385 : 8 : write(sigx1_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol
2386 : 8 : write(sigx1_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha'
2387 : 8 : write(sigx1_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV'
2388 : 8 : write(sigx1_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin'
2389 : 8 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
2390 : 8 : write(sigx1_unt,'(a,i4)') '# Number of core orbitals nc=',nphicor
2391 [ + + ]: 35 : do icor=1,nphicor
2392 [ + + ]: 35 : if (kappacor(icor,itypat_atnbr)==0) then
2393 : 15 : write(sigx1_unt,'(a,2i4,4f15.5)') '# n, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ', &
2394 : 15 : ncor(icor,itypat_atnbr),lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
2395 : 30 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
2396 : : else
2397 [ + + ]: 12 : if (kappacor(icor,itypat_atnbr)>0) then
2398 : 3 : j2=2*lcor(icor,itypat_atnbr)-1
2399 : 3 : write(sigx1_unt,'(a,i4,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
2400 : 3 : & ncor(icor,itypat_atnbr),j2,' /2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
2401 : 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
2402 : : else
2403 [ + + ]: 9 : if (kappacor(icor,itypat_atnbr)<-1) then
2404 : 3 : j2=2*lcor(icor,itypat_atnbr)+1
2405 : 3 : write(sigx1_unt,'(a,i4,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
2406 : 3 : & ncor(icor,itypat_atnbr),j2,'/2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
2407 : 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
2408 : : else
2409 : 6 : write(sigx1_unt,'(a,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
2410 : 6 : & ncor(icor,itypat_atnbr),'1/2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
2411 : 12 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
2412 : : end if
2413 : : end if
2414 : : end if
2415 : : end do
2416 : 8 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
2417 : : write(sigx1_unt,'(a,f10.5,a,f10.5,a)')&
2418 : 8 : & '# Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
2419 : 8 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
2420 [ - + ]: 8 : if(au_units>0) then
2421 : 0 : write(sigx1_unt,'(a)')'# om(au),sig1_av(au),sig1(au) for each orbital'
2422 : : else
2423 : 8 : write(sigx1_unt,'(a)')'# om(eV),sig1_av(Ohm.cm)-1,sig1(Ohm.cm)-1 for each orbital'
2424 : : endif
2425 [ + + ]: 8008 : do iom=1,mom
2426 [ + - ]: 8008 : if(nsppol==1) then
2427 : : write(sigx1_unt,'(100(1x,e15.8))') &
2428 [ + - + + ]: 35000 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1),sigx1(icor,iom,atnbr,1),icor=1,nphicor)
2429 : : else
2430 : : write(sigx1_unt,'(100(1x,e15.8))') &
2431 [ # # # # ]: 0 : & (oml_edge(icor,iom),sum_spin_sigx1_av(icor,iom),sum_spin_sigx1(icor,iom,atnbr),icor=1,nphicor)
2432 : : endif
2433 : : end do
2434 : 8 : close(sigx1_unt)
2435 : : ! absX file
2436 [ + - ]: 8 : if(au_units==0) then
2437 [ - + ]: 8 : if(open_file(trim(filnam_out)//'_absX_at'//str_atm,msg,newunit=absx_unt,form='formatted',action="write")/=0) then
2438 : 0 : ABI_ERROR(msg)
2439 : : end if
2440 : 8 : write(absx_unt,'(a)') '# om(eV), abso X average (cm-1), abso X(cm-1)'
2441 [ + + ]: 8008 : do iom=1,mom
2442 [ + - ]: 8008 : if(nsppol==1) then
2443 : : write(absx_unt,'(100(1x,e15.8))') &
2444 : 27000 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1)/(Speed_Light_SI*eps0),&
2445 [ + - + + ]: 35000 : & sigx1(icor,iom,atnbr,1)/(Speed_Light_SI*eps0),icor=1,nphicor)
2446 : : else
2447 : : write(absx_unt,'(100(1x,e15.8))') &
2448 : 0 : & (oml_edge(icor,iom),sum_spin_sigx1_av(icor,iom)/(Speed_Light_SI*eps0),&
2449 [ # # # # ]: 0 : & sum_spin_sigx1(icor,iom,atnbr)/(Speed_Light_SI*eps0),icor=1,nphicor)
2450 : : endif
2451 : : end do
2452 : 8 : close(absx_unt)
2453 : : endif
2454 : : end if
2455 : :
2456 : : ! _s_sigX
2457 [ - + ]: 8 : if (need_absorption.and.nsppol==2) then
2458 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_sigX_up_at'//str_atm,msg,newunit=sigx1_up_unt,form='formatted',action="write")/=0) then
2459 : 0 : ABI_ERROR(msg)
2460 : : end if
2461 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_sigX_dn_at'//str_atm,msg,newunit=sigx1_dn_unt,form='formatted',action="write")/=0) then
2462 : 0 : ABI_ERROR(msg)
2463 : : end if
2464 [ # # ]: 0 : do iom=1,mom
2465 : : write(sigx1_up_unt,'(100(1x,e15.8))') &
2466 [ # # # # ]: 0 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1),sigx1(icor,iom,atnbr,1),icor=1,nphicor)
2467 : : write(sigx1_dn_unt,'(100(1x,e15.8))') &
2468 [ # # # # ]: 0 : & (oml_edge(icor,iom),sigx1_av(icor,iom,2),sigx1(icor,iom,atnbr,2),icor=1,nphicor)
2469 : : end do
2470 : 0 : close(sigx1_up_unt)
2471 : 0 : close(sigx1_dn_unt)
2472 : : end if
2473 : :
2474 : : ! _emisX file
2475 [ + + ]: 8 : if (need_emissivity) then
2476 [ - + ]: 6 : if (open_file(trim(filnam_out)//'_emisX_at'//str_atm,msg,newunit=ems_unt,form='formatted',action="write")/=0) then
2477 : 0 : ABI_ERROR(msg)
2478 : : end if
2479 [ - + ]: 6 : if(au_units>0) then
2480 : 0 : write(ems_unt,'(a)')'# om(au),sig1_av(au),sig1(au) for each orbital'
2481 : : else
2482 : 6 : write(ems_unt,'(a)')'# om(eV),sig1_av(Ohm.cm)-1,sig1(Ohm.cm)-1 for each orbital'
2483 : : endif
2484 [ + + ]: 6006 : do iom=1,mom
2485 [ + - ]: 6006 : if(nsppol==1) then
2486 : : write(ems_unt,'(3(3(1x,e15.8),2x))') &
2487 [ + - + + ]: 25000 : & (oml_edge(icor,iom),emisx_av(icor,iom,1),emisx(icor,iom,atnbr,1),icor=1,nphicor)
2488 : : else
2489 : : write(ems_unt,'(3(3(1x,e15.8),2x))') &
2490 [ # # # # ]: 0 : & (oml_edge(icor,iom),sum_spin_emisx_av(icor,iom)/dble(natom_atnbr),sum_spin_emisx(icor,iom,atnbr),icor=1,nphicor)
2491 : : endif
2492 : : end do
2493 : 6 : close(ems_unt)
2494 : : end if
2495 : :
2496 : : ! _s_emisX
2497 [ - + ]: 6 : if (need_emissivity.and.nsppol==2) then
2498 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_emisX_up_at'//str_atm,msg,newunit=ems_up_unt,form='formatted',action="write")/=0) then
2499 : 0 : ABI_ERROR(msg)
2500 : : end if
2501 [ # # ]: 0 : if (open_file(trim(filnam_out)//'_emisX_dn_at'//str_atm,msg,newunit=ems_dn_unt,form='formatted',action="write")/=0) then
2502 : 0 : ABI_ERROR(msg)
2503 : : end if
2504 [ # # ]: 0 : do iom=1,mom
2505 : : write(ems_up_unt,'(3(3(1x,e15.8),2x))') &
2506 [ # # # # ]: 0 : & (oml_edge(icor,iom),emisx_av(icor,iom,1)/dble(natom_atnbr),emisx(icor,iom,atnbr,1),icor=1,nphicor)
2507 : : write(ems_dn_unt,'(3(3(1x,e15.8),2x))') &
2508 [ # # # # ]: 0 : (oml_edge(icor,iom),emisx_av(icor,iom,2)/dble(natom_atnbr),emisx(icor,iom,atnbr,2),icor=1,nphicor)
2509 : : end do
2510 : 0 : close(ems_up_unt)
2511 : 0 : close(ems_dn_unt)
2512 : : end if
2513 : :
2514 : : endif ! master node
2515 : :
2516 : : !---------------------------------------------------------------------------------
2517 : : ! End
2518 : :
2519 : : #if defined(HAVE_GPU_MARKERS)
2520 : : call nvtxEndRange()
2521 : : #endif
2522 : : !Release memory space
2523 [ + - ]: 14 : if (need_absorption) then
2524 [ - + ]: 14 : ABI_FREE(sigx1)
2525 [ - + ]: 14 : ABI_FREE(sigx1_av)
2526 [ - + ]: 14 : if (nsppol==2) then
2527 [ # # ]: 0 : ABI_FREE(sum_spin_sigx1)
2528 [ # # ]: 0 : ABI_FREE(sum_spin_sigx1_av)
2529 : : end if
2530 [ - + ]: 14 : ABI_FREE(dom_var1)
2531 [ - + ]: 14 : ABI_FREE(oml_edge)
2532 : : end if
2533 [ + + ]: 14 : if (need_emissivity) then
2534 [ - + ]: 12 : ABI_FREE(emisx)
2535 [ - + ]: 12 : ABI_FREE(emisx_av)
2536 [ - + ]: 12 : if (nsppol==2) then
2537 [ # # ]: 0 : ABI_FREE(sum_spin_emisx)
2538 [ # # ]: 0 : ABI_FREE(sum_spin_emisx_av)
2539 : : end if
2540 [ - + ]: 12 : ABI_FREE(oml_emis)
2541 : : end if
2542 : 14 : ABI_FREE(typat)
2543 [ - + ]: 14 : ABI_FREE(ncor)
2544 [ - + ]: 14 : ABI_FREE(lcor)
2545 [ - + ]: 14 : ABI_FREE(kappacor)
2546 [ - + ]: 14 : ABI_FREE(energy_cor)
2547 [ - + ]: 14 : ABI_FREE(edge)
2548 [ - + ]: 14 : ABI_FREE(eigen0)
2549 [ - + ]: 14 : ABI_FREE(nband)
2550 [ - + ]: 14 : ABI_FREE(occ)
2551 [ - + ]: 14 : ABI_FREE(occ_cor)
2552 [ - + ]: 14 : ABI_FREE(maxocc_cor)
2553 [ - + ]: 14 : ABI_FREE(wtk)
2554 : 14 : call hdr%free()
2555 : 14 : call destroy_mpi_enreg(mpi_enreg)
2556 : :
2557 : : ABI_NVTX_END_RANGE()
2558 : 14 : call timab(1187,2,tsec)
2559 : :
2560 [ - + - + : 28 : end subroutine conducti_paw_core
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - + - +
- + - + -
+ - + - +
- + - + -
+ ]
2561 : : !!***
2562 : :
2563 : : !----------------------------------------------------------------------
2564 : :
2565 : : !!****f* m_conducti/conducti_paw_core_sum
2566 : : !! NAME
2567 : : !! conducti_paw_core_sum
2568 : : !!
2569 : : !! FUNCTION
2570 : : !! Compute the core-valence contributions (X absorption or X emission) of a buffer
2571 : : !! of bands and accumulate them for all core states, frequencies and atoms.
2572 : : !! Parallelized over (atom, frequency, core state): each one sums over all bands.
2573 : : !!
2574 : : !! INPUTS
2575 : : !! nbuf=number of bands in the buffer
2576 : : !! buf_eig(nbuf)=eigenvalue of each band
2577 : : !! buf_wgt(nbuf,nphicor,natom_atnbr)=weight of each band, for each core state and atom
2578 : : !! (pi/ucvol and occupations included, not the 1/omega factor)
2579 : : !! nphicor=number of core states
2580 : : !! natom_atnbr=number of atoms of the selected type
2581 : : !! mom=number of frequencies
2582 : : !! energy_cor(nphicor)=energies of the core states
2583 : : !! oml(nphicor,mom)=frequencies, for each core state
2584 : : !! broad_mode=0 for gaussian broadening, 1 for lorentzian broadening
2585 : : !! dom=broadening width
2586 : : !! [lor_width(nphicor,mom)]=width of the lorentzian broadening (default: dom)
2587 : : !!
2588 : : !! SIDE EFFECTS
2589 : : !! sig_k(nphicor,mom,natom_atnbr)=contributions (accumulated)
2590 : : !!
2591 : : !! SOURCE
2592 : :
2593 [ - + - + ]: 26 : subroutine conducti_paw_core_sum(nbuf,buf_eig,buf_wgt,nphicor,natom_atnbr,mom,energy_cor,oml,&
2594 [ + + + - : 26 : & broad_mode,dom,sig_k,lor_width)
- + ]
2595 : :
2596 : : !Arguments -----------------------------------
2597 : : !scalars
2598 : : integer,intent(in) :: nbuf,nphicor,natom_atnbr,mom,broad_mode
2599 : : real(dp),intent(in) :: dom
2600 : : !arrays
2601 : : real(dp),intent(in) :: buf_eig(:),buf_wgt(:,:,:),energy_cor(nphicor),oml(nphicor,mom)
2602 : : real(dp),intent(in),optional :: lor_width(:,:)
2603 : : real(dp),intent(inout) :: sig_k(nphicor,mom,natom_atnbr)
2604 : :
2605 : : !Local variables-------------------------------
2606 : : !scalars
2607 : : integer :: iatom_atnbr,ib,icor,iom
2608 : : logical :: has_lor_width
2609 : : real(dp) :: acc,arg,arg_max,diff_eig,dirac,dirac_fact,dirac_width,om,width
2610 : :
2611 : : ! *********************************************************************************
2612 : :
2613 : 26 : has_lor_width=present(lor_width)
2614 : 26 : dirac_fact=one/(dom*sqrt(two_pi))
2615 : 26 : dirac_width=one/(sqrt(two)*dom)
2616 : : !Gaussian: dirac=dirac_fact*exp(-arg) is set to zero below 1d-20, i.e. for arg>arg_max
2617 : : ! (skipping exp for these arguments also avoids its slow underflow path)
2618 : 26 : arg_max=log(dirac_fact*1d20)
2619 : :
2620 : : #if defined(HAVE_GPU_MARKERS)
2621 : : call nvtxStartRange("core_sum")
2622 : : #endif
2623 : : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(acc,arg,diff_eig,dirac,ib,om,width)
2624 [ + + ]: 103 : do iatom_atnbr=1,natom_atnbr
2625 [ + + ]: 77103 : do iom=1,mom
2626 [ + + ]: 345077 : do icor=1,nphicor
2627 : 268000 : om=oml(icor,iom)
2628 : 268000 : acc=zero
2629 [ + - ]: 268000 : if (broad_mode==1) then
2630 : 268000 : width=dom
2631 [ + + ]: 268000 : if (has_lor_width) width=lor_width(icor,iom)
2632 [ + + ]: 2412000 : do ib=1,nbuf
2633 : 2144000 : diff_eig=buf_eig(ib)-energy_cor(icor)
2634 : 2144000 : dirac=width/((diff_eig-om)**2+width**2)/pi
2635 [ - + ]: 2144000 : if (dirac<1d-20) dirac=zero
2636 : 2412000 : acc=acc+buf_wgt(ib,icor,iatom_atnbr)*dirac
2637 : : end do
2638 : : else
2639 [ # # ]: 0 : do ib=1,nbuf
2640 : 0 : diff_eig=buf_eig(ib)-energy_cor(icor)
2641 : 0 : arg=((diff_eig-om)*dirac_width)**2
2642 [ # # ]: 0 : if (arg<=arg_max) then
2643 : 0 : dirac=dirac_fact*exp(-arg)
2644 [ # # ]: 0 : if (dirac<1d-20) dirac=zero
2645 : 0 : acc=acc+buf_wgt(ib,icor,iatom_atnbr)*dirac
2646 : : end if
2647 : : end do
2648 : : end if
2649 : 345000 : sig_k(icor,iom,iatom_atnbr)=sig_k(icor,iom,iatom_atnbr)+acc/abs(om)
2650 : : end do
2651 : : end do
2652 : : end do
2653 : : !$OMP END PARALLEL DO
2654 : : #if defined(HAVE_GPU_MARKERS)
2655 : : call nvtxEndRange()
2656 : : #endif
2657 : :
2658 : 26 : end subroutine conducti_paw_core_sum
2659 : : !!***
2660 : :
2661 : : !----------------------------------------------------------------------
2662 : :
2663 : : !!****f* m_conducti/conducti_nc
2664 : : !! NAME
2665 : : !! conducti_nc
2666 : : !!
2667 : : !! FUNCTION
2668 : : !! This program computes the elements of the optical frequency dependent
2669 : : !! conductivity tensor and the conductivity along the three principal axes
2670 : : !! from the Kubo-Greenwood formula.
2671 : : !!
2672 : : !! INPUTS
2673 : : !! (main routine)
2674 : : !!
2675 : : !! OUTPUT
2676 : : !! (main routine)
2677 : : !!
2678 : : !! NOTES
2679 : : !! bantot
2680 : : !! doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy.
2681 : : !! dom=frequency range
2682 : : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
2683 : : !! eigen11(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
2684 : : !! in reciprocal direction 100
2685 : : !! eigen12(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
2686 : : !! in reciprocal direction 010
2687 : : !! eigen13(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
2688 : : !! in reciprocal direction 001
2689 : : !! ecut=kinetic energy planewave cutoff (hartree).
2690 : : !! entropy= entropy associated with the smearing (adimensional)
2691 : : !! fermie= fermi energy (Hartree)
2692 : : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$).
2693 : : !! gmet_inv(3,3)=inverse of reciprocal space metric ($\textrm{bohr}^{2}$).
2694 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
2695 : : !! kin11= Onsager kinetic coeficient=optical conductivity
2696 : : !! kin12= Onsager kinetic coeficient
2697 : : !! kin21= Onsager kinetic coeficient
2698 : : !! kin22= Onsager kinetic coeficient
2699 : : !! Kth=thermal conductivity
2700 : : !! mom=number of frequency for conductivity computation
2701 : : !! mband=maximum number of bands.
2702 : : !! natom = number of atoms in the unit cell.
2703 : : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
2704 : : !! nelect=number of electrons per unit cell
2705 : : !! nkpt=number of k points in the IBZ for this perturbation
2706 : : !! ngfft(3)=integer fft box dimensions.
2707 : : !! nspinor=number of spinorial components of the wavefunctions.
2708 : : !! nsppol=1 for unpolarized, 2 for spin-polarized.
2709 : : !! ntypat = number of atom types.
2710 : : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
2711 : : !! occopt==option for occupancies
2712 : : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
2713 : : !! rprimd(3,3)=real space primitive translations.
2714 : : !! of primitive translations.
2715 : : !! Sth=thermopower
2716 : : !! tsmear=smearing width (or temperature) in Hartree
2717 : : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
2718 : : !! wind=frequency windows for computations of sigma
2719 : : !! wtk(nkpt)=weight assigned to each k point.
2720 : : !! znucl(natom)=atomic number of atoms
2721 : : !! np_sum=noziere-pines sumrule
2722 : : !! cond_kg(mom)=kubo-greenwood conductivity
2723 : : !!
2724 : : !! SOURCE
2725 : :
2726 : :
2727 : 1 : subroutine conducti_nc(filnam,filnam_out)
2728 : :
2729 : : !Arguments -----------------------------------
2730 : : !scalars
2731 : : character(len=fnlen),intent(in) :: filnam,filnam_out
2732 : :
2733 : : !Local variables-------------------------------
2734 : : !scalars
2735 : : integer,parameter :: formeig0=0,formeig1=1
2736 : : integer :: bantot,bd2tot_index,bdtot0_index,bdtot_index
2737 : : integer :: headform,iband,ii,jj,ikpt,iunt
2738 : : integer :: index_1,iom,isppol,jband,l1,l2,mband,mom,natom,nband1
2739 : : integer :: nrot,iomode
2740 : : integer :: nband_k,nkpt,nlign,nrest,nspinor,nsppol,ntypat
2741 : : integer :: occopt,comm
2742 : : integer :: tens_unt,lij_unt,sig_unt,kth_unt,ocond_unt
2743 : : real(dp) :: deltae,dosdeltae,diff_occ,dom,ecut,entropy,fermie,maxocc
2744 : : real(dp) :: nelect,np_sum,np_sum_k1,np_sum_k2,omin,oml,socc,socc_k,sig
2745 : : real(dp) :: tphysel,tsmear,ucvol,wind,Tatm
2746 : : character(len=fnlen) :: filnam0,filnam1,filnam2,filnam3
2747 : : character(len=500) :: msg
2748 : 1 : type(hdr_type) :: hdr
2749 [ + + + + : 14 : type(wfk_t) :: gswfk,ddk1,ddk2,ddk3
+ + + + ]
2750 : : !arrays
2751 : 1 : integer,allocatable :: nband(:)
2752 : : real(dp) :: gmet(3,3),gmet_inv(3,3),gprimd(3,3),gprimd_inv(3,3),rmet(3,3),rprimd(3,3)
2753 : 2 : real(dp),allocatable :: cond_kg(:,:,:),cond_kg_cart(:,:,:),cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:)
2754 : 3 : real(dp),allocatable :: doccde(:),doccde_k(:),cond_kg_xx(:),cond_kg_yy(:),cond_kg_zz(:),trace(:)
2755 : 1 : real(dp),allocatable :: eig0_k(:),eig0tmp(:),eig1_k(:,:),eigen0(:),eigen11(:)
2756 : 1 : real(dp),allocatable :: eigen12(:),eigtmp(:)
2757 : 2 : real(dp),allocatable :: eigen13(:),occ(:),occ_k(:),wtk(:),cond_tot(:),oml1(:)
2758 : 1 : real(dp),allocatable :: kin11(:),kin12(:),kin21(:),kin22(:)
2759 : 1 : real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
2760 : : real(dp) :: cond_kg_w(3,3),z(3,3)
2761 : : real(dp) :: eig_cond(3)
2762 : :
2763 : : ! *********************************************************************************
2764 : :
2765 : : ABI_NVTX_START_RANGE(NVTX_CONDUCTI_NC)
2766 : :
2767 : : !Read data file
2768 [ - + ]: 1 : if (open_file(filnam,msg,newunit=iunt,form='formatted',status="old")/=0) then
2769 : 0 : ABI_ERROR(msg)
2770 : : end if
2771 : :
2772 : 1 : rewind(iunt)
2773 : 1 : read(iunt,*)
2774 : 1 : read(iunt,'(a)')filnam1 ! first ddk file
2775 : 1 : read(iunt,'(a)')filnam2 ! second ddk file
2776 : 1 : read(iunt,'(a)')filnam3 ! third ddk file
2777 : 1 : read(iunt,'(a)')filnam0 ! ground-state data
2778 : :
2779 : : !Open the GS Wavefunction file and the 3 DDK files.
2780 : : ! TODO: one should perform basic consistency tests for the GS WFK and the DDK files, e.g.
2781 : : ! k-points and their order, spins, number of bands could differ in the four files.
2782 : : ! Note indeed that we are assuming the same numer of bands in all the files.
2783 : 1 : comm = xmpi_comm_self
2784 : 1 : call nctk_fort_or_ncfile(filnam0, iomode, msg)
2785 [ - + ]: 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
2786 : 1 : call gswfk%open_read(filnam0, formeig0, iomode, get_unit(), comm)
2787 : :
2788 : 1 : call nctk_fort_or_ncfile(filnam1, iomode, msg)
2789 [ - + ]: 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
2790 : 1 : call ddk1%open_read(filnam1, formeig1, iomode, get_unit(), comm, hdr_out=hdr)
2791 : :
2792 : 1 : call nctk_fort_or_ncfile(filnam2, iomode, msg)
2793 [ - + ]: 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
2794 : 1 : call ddk2%open_read(filnam2, formeig1, iomode, get_unit(), comm)
2795 : :
2796 : 1 : call nctk_fort_or_ncfile(filnam3, iomode, msg)
2797 [ - + ]: 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
2798 : 1 : call ddk3%open_read(filnam3, formeig1, iomode, get_unit(), comm)
2799 : :
2800 [ - + ]: 1 : if (ddk1%compare(ddk2) /= 0) then
2801 : 0 : ABI_ERROR("ddk1 and ddk2 are not consistent. see above messages")
2802 : : end if
2803 [ - + ]: 1 : if (ddk1%compare(ddk3) /= 0) then
2804 : 0 : ABI_ERROR("ddk1 and ddk3 are not consistent. see above messages")
2805 : : end if
2806 : :
2807 : : !Extract params from the header of the first ddk file (might have been the GS file ?)
2808 : :
2809 : : !Extract info from the header
2810 : 1 : headform=hdr%headform
2811 : 1 : bantot=hdr%bantot
2812 : 1 : ecut=hdr%ecut_eff
2813 : 1 : natom=hdr%natom
2814 : 1 : nkpt=hdr%nkpt
2815 : 1 : nspinor=hdr%nspinor
2816 : 1 : nsppol=hdr%nsppol
2817 : 1 : ntypat=hdr%ntypat
2818 : 1 : occopt=hdr%occopt
2819 [ + + + + ]: 13 : rprimd(:,:)=hdr%rprimd(:,:)
2820 [ - + + - : 3 : ABI_MALLOC(nband,(nkpt*nsppol))
+ - + - -
+ - + ]
2821 [ - + + - : 3 : ABI_MALLOC(occ,(bantot))
+ - + - -
+ - + ]
2822 : 1 : fermie=hdr%fermie
2823 [ + + ]: 21 : occ(1:bantot)=hdr%occ(1:bantot)
2824 [ + + ]: 2 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
2825 : :
2826 : : !Get mband, as the maximum value of nband(nkpt)
2827 [ + + ]: 2 : mband=maxval(nband(:))
2828 : :
2829 : 1 : write(std_out,*)
2830 : 1 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
2831 : 1 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
2832 : 1 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
2833 : 1 : write(std_out,'(a,i8)') ' natom =',natom
2834 : 1 : write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband
2835 : 1 : write(std_out,'(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
2836 : 1 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
2837 : :
2838 : : !Prepare the reading of ddk Wff files
2839 [ + - - + : 3 : ABI_MALLOC(eigtmp,(2*mband*mband))
+ - + - +
- - + -
+ ]
2840 [ - + + - : 3 : ABI_MALLOC(eig0tmp,(mband))
+ - + - -
+ - + ]
2841 : :
2842 : : !Read the eigenvalues of ground-state and ddk files
2843 [ - + + - : 3 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
+ - + - -
+ - + ]
2844 [ + - - + : 3 : ABI_MALLOC(eigen11,(2*mband*mband*nkpt*nsppol))
+ - + - +
- - + -
+ ]
2845 [ - + + - : 2 : ABI_MALLOC(eigen12,(2*mband*mband*nkpt*nsppol))
- + - + ]
2846 [ - + + - : 2 : ABI_MALLOC(eigen13,(2*mband*mband*nkpt*nsppol))
- + - + ]
2847 : 1 : bdtot0_index=0 ; bdtot_index=0
2848 [ + + ]: 2 : do isppol=1,nsppol
2849 [ + + ]: 3 : do ikpt=1,nkpt
2850 : 1 : nband1=nband(ikpt+(isppol-1)*nkpt)
2851 : 1 : call gswfk%read_eigk(ikpt,isppol,xmpio_single,eig0tmp)
2852 [ + + ]: 21 : eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1)
2853 : :
2854 : 1 : call ddk1%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2855 [ + + ]: 801 : eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2856 : :
2857 : 1 : call ddk2%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2858 [ + + ]: 801 : eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2859 : :
2860 : 1 : call ddk3%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2861 [ + + ]: 801 : eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2862 : :
2863 : 1 : bdtot0_index=bdtot0_index+nband1
2864 : 2 : bdtot_index=bdtot_index+2*nband1**2
2865 : : end do
2866 : : end do
2867 : :
2868 : : !Close files
2869 : 1 : call gswfk%close()
2870 : 1 : call ddk1%close()
2871 : 1 : call ddk2%close()
2872 : 1 : call ddk3%close()
2873 : :
2874 [ - + ]: 1 : ABI_FREE(eigtmp)
2875 [ - + ]: 1 : ABI_FREE(eig0tmp)
2876 : :
2877 : : !---------------------------------------------------------------------------------
2878 : : !Gmet inversion
2879 : 1 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2880 : 1 : call matr3inv(gmet,gmet_inv)
2881 : 1 : call matr3inv(gprimd,gprimd_inv)
2882 : :
2883 : : !---------------------------------------------------------------------------------
2884 : : !Derivative of occupation wrt the energy.
2885 : :
2886 [ - + + - : 2 : ABI_MALLOC(doccde,(mband*nkpt*nsppol))
- + - + ]
2887 [ - + + - : 3 : ABI_MALLOC(wtk,(nkpt))
+ - + - -
+ - + ]
2888 : :
2889 : 1 : read(iunt,*)tsmear
2890 : 1 : Tatm=tsmear*Ha_K
2891 : 1 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
2892 : : !
2893 : 1 : nlign=nkpt/6
2894 : 1 : nrest=nkpt-6*nlign
2895 : 1 : index_1=0
2896 [ - + ]: 1 : do ii=1,nlign
2897 : 0 : read(iunt,*)wtk(1+index_1:6+index_1)
2898 : 1 : index_1=index_1+6
2899 : : end do
2900 [ + - ]: 1 : if (nrest/=0) then
2901 : 1 : read(iunt,*)wtk(6*nlign+1:nkpt)
2902 : : end if
2903 : : !
2904 [ - + ]: 1 : if (occopt==1) then
2905 : 0 : write(std_out,'(a,i4)') ' occopt =',occopt
2906 [ # # ]: 0 : doccde=zero
2907 : : else
2908 : 1 : tphysel=zero
2909 : 1 : maxocc=two/(nsppol*nspinor)
2910 : 1 : dosdeltae=zero
2911 : : ! CP: using 1 and nband(0) as dummy value, because function
2912 : : ! not implemented for occopt==9; adding fermih=fermie in the list of arguments as well
2913 : : call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermie,maxocc,mband,nband,&
2914 : 1 : & nelect,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,11,wtk,1,nband(1))
2915 : : end if
2916 : :
2917 : : !---------------------------------------------------------------------------------
2918 : : !Size of the frequency range
2919 : :
2920 : 1 : read(iunt,*)dom,wind
2921 : 1 : close(iunt)
2922 : 1 : mom=int(wind/dom)
2923 [ - + + - : 3 : ABI_MALLOC(oml1,(mom))
+ - + - -
+ - + ]
2924 [ + + ]: 273 : do iom=1,mom
2925 : 273 : oml1(iom)=tol10*1000._dp+dble(iom)*dom
2926 : : end do
2927 : :
2928 [ - + + - : 4 : ABI_MALLOC(cond_nd,(mom,3,3))
+ - + - +
- - + -
+ ]
2929 [ - + + - : 2 : ABI_MALLOC(cond_kg,(mom,3,3))
- + - + ]
2930 [ - + + - : 2 : ABI_MALLOC(cond_kg_cart,(mom,3,3))
- + - + ]
2931 [ - + + - : 2 : ABI_MALLOC(cond_kg_xx,(mom))
- + - + ]
2932 [ - + + - : 2 : ABI_MALLOC(cond_kg_yy,(mom))
- + - + ]
2933 [ - + + - : 2 : ABI_MALLOC(trace,(mom))
- + - + ]
2934 [ - + + - : 2 : ABI_MALLOC(cond_kg_zz,(mom))
- + - + ]
2935 [ - + + - : 2 : ABI_MALLOC(cond_tot,(mom))
- + - + ]
2936 [ - + + - : 2 : ABI_MALLOC(kin11,(mom))
- + - + ]
2937 [ - + + - : 2 : ABI_MALLOC(kin12,(mom))
- + - + ]
2938 [ - + + - : 2 : ABI_MALLOC(kin21,(mom))
- + - + ]
2939 [ - + + - : 2 : ABI_MALLOC(kin22,(mom))
- + - + ]
2940 [ - + + - : 2 : ABI_MALLOC(kin11_k,(mom))
- + - + ]
2941 [ - + + - : 2 : ABI_MALLOC(kin12_k,(mom))
- + - + ]
2942 [ - + + - : 2 : ABI_MALLOC(kin21_k,(mom))
- + - + ]
2943 [ - + + - : 2 : ABI_MALLOC(kin22_k,(mom))
- + - + ]
2944 [ - + + - : 2 : ABI_MALLOC(Kth,(mom))
- + - + ]
2945 [ - + + - : 2 : ABI_MALLOC(Stp,(mom))
- + - + ]
2946 : 1 : write(std_out,'(a,i8,2f10.5,a)')' mom,wind,dom =',mom,wind,dom,' Ha'
2947 : :
2948 : : !---------------------------------------------------------------------------------
2949 : :
2950 [ + + ]: 273 : kin11 = zero
2951 [ + + ]: 273 : kin12 = zero
2952 [ + + ]: 273 : kin21 = zero
2953 [ + + ]: 273 : kin22 = zero
2954 : 1 : np_sum = zero
2955 : 1 : socc = zero
2956 [ + + + + : 2461 : cond_kg = zero
+ + ]
2957 : :
2958 : : !LOOP OVER SPINS
2959 [ + + ]: 2 : do isppol=1,nsppol
2960 : : !
2961 : 1 : bdtot_index = 0
2962 : 1 : bd2tot_index = 0
2963 : :
2964 : 1 : deltae = zero
2965 : : !
2966 : : ! BIG FAT k POINT LOOP
2967 : : !
2968 [ + + ]: 2 : do ikpt=1,nkpt
2969 : :
2970 : 1 : nband_k=nband(ikpt+(isppol-1)*nkpt)
2971 : :
2972 [ - + + - : 3 : ABI_MALLOC(eig0_k,(nband_k))
+ - + - -
+ - + ]
2973 [ + - - + : 4 : ABI_MALLOC(eig1_k,(2*nband_k**2,3))
+ - + - +
- + - - +
- + ]
2974 [ - + + - : 2 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
2975 [ - + + - : 2 : ABI_MALLOC(doccde_k,(nband_k))
- + - + ]
2976 [ - + + - : 4 : ABI_MALLOC(dhdk2_r,(3,3,nband_k,nband_k))
- + + - +
- + - - +
- + ]
2977 [ - + + - : 4 : ABI_MALLOC(dhdk2_g,(nband_k,nband_k))
- + + - +
- + - - +
- + ]
2978 : :
2979 [ + + + + : 2461 : cond_nd = zero
+ + ]
2980 [ + + ]: 273 : kin11_k = zero
2981 [ + + ]: 273 : kin12_k = zero
2982 [ + + ]: 273 : kin21_k = zero
2983 [ + + ]: 273 : kin22_k = zero
2984 : 21 : np_sum_k1 = zero
2985 : 21 : np_sum_k2 = zero
2986 : 21 : socc_k = zero
2987 [ + + + + : 5221 : dhdk2_r = zero
+ + + + ]
2988 [ + + + + ]: 421 : dhdk2_g = zero
2989 : :
2990 : : ! eigenvalue for k-point
2991 [ + + ]: 21 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
2992 : : ! first derivative eigenvalues for k-point
2993 [ + + ]: 801 : eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2994 [ + + ]: 801 : eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2995 [ + + ]: 801 : eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2996 : : ! occupation numbers for k-point
2997 [ + + ]: 21 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
2998 : : ! derivative of occupation number for k-point
2999 [ + + ]: 21 : doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
3000 : :
3001 : : ! LOOP OVER BAND
3002 [ + + ]: 21 : do iband=1,nband_k
3003 [ + + ]: 420 : do jband=1,nband_k
3004 : : !
3005 : : ! TODO : replace with BLAS calls
3006 [ + + ]: 1600 : do l1=1,3
3007 [ + + ]: 5200 : do l2=1,3
3008 [ + + ]: 15600 : do ii=1,3
3009 [ + + ]: 46800 : do jj=1,3
3010 : : dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)+(rprimd(l1,ii)&
3011 : : & *eig1_k(2*iband-1+(jband-1)*2*nband_k,ii)*&
3012 : : & rprimd(l2,jj)*eig1_k(2*iband-1+(jband-1)*2*nband_k,jj)&
3013 : : & +rprimd(l1,ii)*eig1_k(2*iband +(jband-1)*2*nband_k,ii)*&
3014 : 43200 : & rprimd(l2,jj)*eig1_k(2*iband+(jband-1)*2*nband_k,jj))
3015 : : end do
3016 : : end do
3017 : : end do
3018 : : end do
3019 : :
3020 [ + + ]: 1600 : do l1=1,3
3021 [ + + ]: 5200 : do l2=1,3
3022 : 4800 : dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)/two_pi/two_pi
3023 : : end do
3024 : : end do
3025 : :
3026 : : ! TODO: replace with BLAS calls
3027 [ + + ]: 1600 : do l1=1,3
3028 [ + + ]: 5200 : do l2=1,3
3029 : : dhdk2_g(iband,jband)=dhdk2_g(iband,jband)+gmet_inv(l1,l2)*( &
3030 : : & eig1_k(2*iband-1+(jband-1)*2*nband_k,l1)*&
3031 : : & eig1_k(2*iband-1+(jband-1)*2*nband_k,l2) &
3032 : : & +eig1_k(2*iband +(jband-1)*2*nband_k,l1)*&
3033 : 4800 : & eig1_k(2*iband +(jband-1)*2*nband_k,l2))
3034 : : end do
3035 : : end do
3036 : 400 : dhdk2_g(iband,jband)=dhdk2_g(iband,jband)/two_pi/two_pi
3037 : :
3038 : 400 : diff_occ = occ_k(iband)-occ_k(jband)
3039 : : ! if (dabs(diff_occ)>=tol8) then
3040 : :
3041 : : ! Conductivity for each omega
3042 : 400 : omin = zero
3043 [ + + ]: 109200 : do iom=1,mom
3044 : 108800 : oml=oml1(iom)
3045 [ + + ]: 108800 : if (jband>iband) then
3046 : : sig= dhdk2_g(iband,jband)&
3047 : : & *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)&
3048 : 51680 : & -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2))
3049 : 51680 : kin11_k(iom)=kin11_k(iom)+sig
3050 : 51680 : kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie)
3051 : 51680 : kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie)
3052 : : kin22_k(iom)=kin22_k(iom) + &
3053 : 51680 : & sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie)
3054 : : end if
3055 [ + + ]: 435600 : do l1=1,3
3056 [ + + ]: 1414400 : do l2=1,3
3057 : : cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(l1,l2,iband,jband)&
3058 : 1305600 : & *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)
3059 : : end do
3060 : : end do
3061 : :
3062 : : end do
3063 : :
3064 : : ! Sumrule start
3065 [ + + ]: 420 : if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then
3066 : : np_sum_k1=np_sum_k1 -dhdk2_g(iband,jband)&
3067 : 380 : & *(diff_occ)/(eig0_k(iband)-eig0_k(jband))
3068 : : else
3069 : 20 : np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(iband,jband)
3070 : : end if
3071 : :
3072 : :
3073 : : ! end loop over band
3074 : : end do
3075 : 21 : socc_k=socc_k+occ_k(iband)
3076 : : end do
3077 : :
3078 [ + + ]: 273 : do iom=1,mom
3079 : 272 : kin11(iom)=kin11(iom)+wtk(ikpt)*kin11_k(iom)
3080 : 272 : kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
3081 : 272 : kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
3082 : 272 : kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
3083 [ + + ]: 1089 : do l1=1,3
3084 [ + + ]: 3536 : do l2=1,3
3085 : 3264 : cond_kg(iom,l1,l2)=cond_kg(iom,l1,l2)+wtk(ikpt)*cond_nd(iom,l1,l2)
3086 : : end do
3087 : : end do
3088 : : end do
3089 : :
3090 : 1 : np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2)
3091 : 1 : socc=socc+wtk(ikpt)*socc_k
3092 : :
3093 : : ! Validity limit
3094 : 1 : deltae=deltae+(eig0_k(nband_k)-fermie)
3095 : :
3096 : 1 : bd2tot_index=bd2tot_index+2*nband_k**2
3097 : 1 : bdtot_index=bdtot_index+nband_k
3098 : 1 : ABI_FREE(eig0_k)
3099 [ - + ]: 1 : ABI_FREE(eig1_k)
3100 [ - + ]: 1 : ABI_FREE(occ_k)
3101 [ - + ]: 1 : ABI_FREE(doccde_k)
3102 [ - + ]: 1 : ABI_FREE(dhdk2_r)
3103 [ - + ]: 2 : ABI_FREE(dhdk2_g)
3104 : : ! End loop over k
3105 : : end do
3106 : :
3107 : 1 : write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three,socc
3108 : : write(std_out,'(a,f10.5,a,f10.5,a)')&
3109 : 2 : & ' Emax-Efermi =',deltae/dble(nkpt),' Ha',deltae/dble(nkpt)*Ha_eV,' eV'
3110 : :
3111 : : !End loop over spins
3112 : : end do
3113 : :
3114 [ + + + + : 2461 : cond_kg=cond_kg*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
+ + ]
3115 : :
3116 : :
3117 : : !Check that new output file does NOT exist
3118 : : !Keep this line : prevent silly (compiler ?) bug on HP 8000
3119 : 1 : write(std_out,*)' conducti : call isfile '
3120 : : !
3121 [ - + ]: 1 : if (open_file(trim(filnam_out)//'_tens',msg,newunit=tens_unt,form='formatted',action="write")/=0) then
3122 : 0 : ABI_ERROR(msg)
3123 : : end if
3124 [ - + ]: 1 : if (open_file(trim(filnam_out)//'_Lij',msg,newunit=lij_unt,form='formatted',action="write")/=0) then
3125 : 0 : ABI_ERROR(msg)
3126 : : end if
3127 : 1 : write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22'
3128 : :
3129 [ - + ]: 1 : if (open_file(trim(filnam_out)//'_sig',msg,newunit=sig_unt,form='formatted',action="write")/=0) then
3130 : 0 : ABI_ERROR(msg)
3131 : : end if
3132 : 1 : write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1'
3133 : :
3134 [ - + ]: 1 : if (open_file(trim(filnam_out)//'_Kth',msg,newunit=kth_unt,form='formatted',action="write")/=0) then
3135 : 0 : ABI_ERROR(msg)
3136 : : end if
3137 : : write(kth_unt,'(a)')&
3138 : 1 : & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)'
3139 : :
3140 [ - + ]: 1 : if (open_file(trim(filnam_out)//'.out',msg,newunit=ocond_unt,form='formatted',action="write")/=0) then
3141 : 0 : ABI_ERROR(msg)
3142 : : end if
3143 : 1 : write(ocond_unt,'(a)' )' Conducti output file:'
3144 : 1 : write(ocond_unt,'(a)' )' Contains all results produced by conducti utility'
3145 : 1 : write(ocond_unt,'(a)' )' '
3146 : 1 : write(ocond_unt,'(a)')' # omega(ua) cond(ua) thermal cond(ua) thermopower(ua)'
3147 : :
3148 : : !Keep this line : prevent silly (compiler ?) bug on HP 8000
3149 : 1 : write(std_out,*)' conducti : after call isfile '
3150 : :
3151 : : !Compute thermal conductivity and thermopower
3152 [ + + ]: 273 : do iom=1,mom
3153 : 272 : oml=oml1(iom)
3154 : 272 : kin11(iom)=kin11(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
3155 : 272 : kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
3156 : 272 : kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
3157 : 272 : kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
3158 [ + + ]: 272 : if (dabs(kin11(iom))<10.0d-20) kin11(iom)=zero
3159 : 272 : Kth(iom)=kin22(iom)
3160 : 272 : Stp(iom)=zero
3161 [ + + ]: 272 : if(kin11(iom)/=zero) then
3162 : 107 : Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/kin11(iom))
3163 : 107 : Stp(iom)=kin12(iom)/(kin11(iom)*Tatm)
3164 : : end if
3165 [ + + ]: 272 : if (dabs(Kth(iom))<10.0d-20) Kth(iom)=zero
3166 [ + + ]: 272 : if (dabs(Stp(iom))<10.0d-20) Stp(iom)=zero
3167 [ + + ]: 272 : if (dabs(kin12(iom))<10.0d-20) kin12(iom)=zero
3168 [ + + ]: 272 : if (dabs(kin21(iom))<10.0d-20) kin21(iom)=zero
3169 [ + + ]: 272 : if (dabs(kin22(iom))<10.0d-20) kin22(iom)=zero
3170 : :
3171 : 272 : write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9
3172 : 272 : write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,kin11(iom),kin11(iom)*Ohmcm
3173 : 272 : write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm,Stp(iom),Stp(iom)*3.6753d-2
3174 : 273 : write(ocond_unt,'(1f12.5,3es22.12)') oml,kin11(iom),Kth(iom),Stp(iom)
3175 : : end do
3176 : :
3177 : 1 : write(tens_unt,'(a)' )' Conductivity file '
3178 : 1 : write(tens_unt,'(a)' )' ----------------- '
3179 : 1 : write(tens_unt,'(a)' )' Contain first the full conductivity tensor, for the desired set of energies,'
3180 : 1 : write(tens_unt,'(a)' )' then, the three principal values, for the desired set of energies'
3181 : 1 : write(tens_unt,'(a)' )' (note that eigenvalues are not directly associated with xx,yy,zz)'
3182 : 1 : write(tens_unt,'(a)' )' '
3183 : :
3184 : 1 : write(ocond_unt,'(a)' )' '
3185 : 1 : write(ocond_unt,'(a)' )' full conductivity tensor, for the desired set of energies'
3186 : 1 : write(ocond_unt,'(a)' )' then, the three principal values, for the desired set of energies:'
3187 : :
3188 [ + + ]: 273 : do iom=1,mom
3189 : 272 : oml=oml1(iom)*Ha_eV
3190 : 272 : write(tens_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
3191 : 272 : write(ocond_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
3192 [ + + ]: 1089 : do l1=1,3
3193 [ + - + + ]: 3264 : write(tens_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
3194 [ + - + + ]: 3536 : write(ocond_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
3195 : : end do
3196 : : end do
3197 : :
3198 : : !Diagonalizing the conductivity matrix for sigma_xx,sigma_yy,sigma_zz
3199 [ + + ]: 273 : cond_kg_xx=0d0
3200 [ + + ]: 273 : cond_kg_yy=0d0
3201 [ + + ]: 273 : cond_kg_zz=0d0
3202 : : !trace=0d0 ! Used for checking with the original version of the code
3203 [ + + ]: 273 : do iom=1,mom
3204 : 272 : oml=oml1(iom)*Ha_eV
3205 : 272 : cond_kg_w=0d0
3206 [ + + ]: 1088 : do l1=1,3
3207 [ + + ]: 3536 : do l2=1,3
3208 : 3264 : cond_kg_w(l1,l2)=cond_kg(iom,l1,l2)
3209 : : end do
3210 : : end do
3211 : 272 : call jacobi(cond_kg_w,3,3,eig_cond,z,nrot)
3212 : :
3213 : : ! When the value is too small, set it to zero before printing
3214 [ + + ]: 272 : if(abs(eig_cond(1))<tol10)eig_cond(1)=zero
3215 [ + + ]: 272 : if(abs(eig_cond(2))<tol10)eig_cond(2)=zero
3216 [ + + ]: 272 : if(abs(eig_cond(3))<tol10)eig_cond(3)=zero
3217 : :
3218 : 272 : cond_kg_xx(iom)=eig_cond(1)
3219 : 272 : cond_kg_yy(iom)=eig_cond(2)
3220 : 273 : cond_kg_zz(iom)=eig_cond(3)
3221 : : ! trace(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
3222 : : end do
3223 : :
3224 : : !DEBUG Keep this line : prevent silly (compiler ?) bug on HP 8000
3225 : : !write(std_out,*)' conducti : after open '
3226 : : !ENDDEBUG
3227 : :
3228 : 1 : write(tens_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
3229 : 1 : write(tens_unt,'(a)')' '
3230 : 1 : write(tens_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)'
3231 : :
3232 : 1 : write(ocond_unt,'(a)')' '
3233 : 1 : write(ocond_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
3234 : 1 : write(ocond_unt,'(a)')' '
3235 : 1 : write(ocond_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)'
3236 : :
3237 : :
3238 [ + + ]: 273 : do iom=1,mom
3239 : 272 : cond_tot(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
3240 : 272 : write(tens_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
3241 : 273 : write(ocond_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
3242 : : end do
3243 : :
3244 : 1 : write(tens_unt,*)
3245 : 1 : write(tens_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1'
3246 : 1 : write(ocond_unt,*)
3247 : 1 : write(ocond_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1'
3248 : :
3249 [ + + ]: 273 : do iom=1,mom
3250 : 272 : oml=oml1(iom)*Ha_eV
3251 : 272 : cond_tot(iom)=cond_tot(iom)*Ohmcm
3252 : 272 : cond_kg_xx(iom)=cond_kg_xx(iom)*Ohmcm
3253 : 272 : cond_kg_yy(iom)=cond_kg_yy(iom)*Ohmcm
3254 : 272 : cond_kg_zz(iom)=cond_kg_zz(iom)*Ohmcm
3255 : 272 : write(tens_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
3256 : 273 : write(ocond_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
3257 : : end do
3258 : :
3259 : : !Calculate the imaginary part of the conductivity (principal value)
3260 : : !+derived optical properties.
3261 : 1 : call msig(kin11,mom,oml1,filnam_out,zero,0)
3262 : :
3263 : 1 : close(tens_unt)
3264 : 1 : close(lij_unt)
3265 : 1 : close(sig_unt)
3266 : 1 : close(kth_unt)
3267 : 1 : close(ocond_unt)
3268 : :
3269 [ - + ]: 1 : ABI_FREE(nband)
3270 [ - + ]: 1 : ABI_FREE(oml1)
3271 [ - + ]: 1 : ABI_FREE(occ)
3272 [ - + ]: 1 : ABI_FREE(eigen11)
3273 [ - + ]: 1 : ABI_FREE(eigen12)
3274 [ - + ]: 1 : ABI_FREE(eigen13)
3275 [ - + ]: 1 : ABI_FREE(eigen0)
3276 [ - + ]: 1 : ABI_FREE(doccde)
3277 [ - + ]: 1 : ABI_FREE(wtk)
3278 [ - + ]: 1 : ABI_FREE(cond_nd)
3279 [ - + ]: 1 : ABI_FREE(cond_kg)
3280 [ - + ]: 1 : ABI_FREE(cond_kg_cart)
3281 [ - + ]: 1 : ABI_FREE(cond_kg_xx)
3282 [ - + ]: 1 : ABI_FREE(cond_kg_yy)
3283 [ - + ]: 1 : ABI_FREE(trace)
3284 [ - + ]: 1 : ABI_FREE(cond_kg_zz)
3285 [ - + ]: 1 : ABI_FREE(cond_tot)
3286 [ - + ]: 1 : ABI_FREE(kin11)
3287 [ - + ]: 1 : ABI_FREE(kin22)
3288 [ - + ]: 1 : ABI_FREE(kin12)
3289 [ - + ]: 1 : ABI_FREE(kin21)
3290 [ - + ]: 1 : ABI_FREE(kin11_k)
3291 [ - + ]: 1 : ABI_FREE(kin22_k)
3292 [ - + ]: 1 : ABI_FREE(kin12_k)
3293 [ - + ]: 1 : ABI_FREE(kin21_k)
3294 [ - + ]: 1 : ABI_FREE(Stp)
3295 [ - + ]: 1 : ABI_FREE(Kth)
3296 : 1 : call hdr%free()
3297 : :
3298 : : ABI_NVTX_END_RANGE()
3299 : :
3300 [ - + - + : 5 : end subroutine conducti_nc
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - - - -
- - - - -
+ - + -
+ ]
3301 : : !!***
3302 : :
3303 : : !----------------------------------------------------------------------
3304 : :
3305 : : !!****f* m_conducti/msig
3306 : : !! NAME
3307 : : !! msig
3308 : : !!
3309 : : !! FUNCTION
3310 : : !! This program computes the elements of the optical frequency dependent
3311 : : !! conductivity tensor and the conductivity along the three principal axes
3312 : : !! from the Kubo-Greenwood formula for PAW formalism
3313 : : !!
3314 : : !! INPUTS
3315 : : !! fcti(npti)= conductivity, as calculated in conducti
3316 : : !! npti= number of points to calculate conductivity
3317 : : !! xi(npti)= energies where the conductivity is calculated
3318 : : !!
3319 : : !! OUTPUT
3320 : : !! no output, only files
3321 : : !!
3322 : : !! NOTES
3323 : : !! this program calculates the imaginary part of the conductivity (principal value)
3324 : : !! +derived optical properties.
3325 : : !!
3326 : : !! SOURCE
3327 : :
3328 : 9 : subroutine msig(fcti,npti,xi,filnam_out_sig,phi,au_units)
3329 : :
3330 : : !Arguments -----------------------------------
3331 : : !scalars
3332 : : integer,intent(in) :: npti, au_units
3333 : : !arrays
3334 : : real(dp),intent(in) :: fcti(npti),xi(npti)
3335 : : character(len=fnlen),intent(in) :: filnam_out_sig
3336 : : real(dp),intent(in) :: phi
3337 : :
3338 : : !Local variables-------------------------------
3339 : : !scalars
3340 : : integer :: ii,ip,eps_unt,abs_unt
3341 : : real(dp),parameter :: del=0.001_dp
3342 : : real(dp) :: dx,eps1,eps2,komega,pole,refl_s,refl_p,sigma2,xsum,ff,ffp,ffpp,abso,sigma1
3343 : : character(len=500) :: msg
3344 : : !arrays
3345 : 9 : real(dp),allocatable :: fct(:)
3346 : : real(dp) :: xx1,xx2,xx,nomega
3347 : : complex(dp) :: epsc,cos_phi,sin_phi,sqroot,crefl_s,crefl_p
3348 : : ! *********************************************************************************
3349 : :
3350 : :
3351 : 9 : write(std_out,'(2a)')ch10,'Calculate the principal value and related optical properties'
3352 : 9 : write(std_out,'(a)')'Use default value for delta interval: del=1e-3'
3353 : :
3354 [ - + ]: 9 : if (open_file(trim(filnam_out_sig)//'_eps',msg,newunit=eps_unt,status='replace',action="write")/=0) then
3355 : 0 : ABI_ERROR(msg)
3356 : : end if
3357 [ + - ]: 9 : if(au_units==0) then
3358 : 9 : write(eps_unt,'(a)')'#energy (eV),sigma_1(Ohm-1cm-1),sigma_2(Ohm-1cm-1),epsilon_1(cgs),epsilon_2(cgs)'
3359 : : else
3360 : 0 : write(eps_unt,'(a)')'#energy (Ha),sigma_1(au),sigma_2(au),epsilon_1(au),epsilon_2(au)'
3361 : : endif
3362 : :
3363 [ - + ]: 9 : if (open_file(trim(filnam_out_sig)//'_abs',msg,newunit=abs_unt,status='replace',action="write")/=0) then
3364 : 0 : ABI_ERROR(msg)
3365 : : end if
3366 [ + - ]: 9 : if(au_units==0) then
3367 : 9 : write(abs_unt,'(a)')'#energy(eV),nomega,komega,refl. s,refl. p,abso.(cm-1)'
3368 : : else
3369 : 0 : write(abs_unt,'(a)')'#energy(Ha),nomega,komega,refl. s,refl. p,abso.(au)'
3370 : : endif
3371 : :
3372 [ - + + - : 27 : ABI_MALLOC(fct,(npti))
+ - + - -
+ - + ]
3373 : :
3374 : : !loop on the initial energy grid
3375 [ + + ]: 8281 : do ip=1,npti
3376 : : !!! Taylor expansion up to second order of fcti at xx=xi(ii) for each ii
3377 : : !!! Then, in each discretization interval, the integral can be performed analytically
3378 : 8272 : pole=xi(ip)
3379 : 8272 : dx=(xi(npti)-xi(1))/dble(npti-1)
3380 : 8272 : xsum=zero
3381 [ + + ]: 8082256 : do ii=1,npti
3382 : 8073984 : xx=xi(ii)
3383 : 8082256 : fct(ii)=fcti(ii)*pole/(xx+pole)
3384 : : enddo
3385 [ + + ]: 8082256 : do ii=1,npti
3386 : 8073984 : xx=xi(ii)
3387 : 8073984 : ff=fct(ii)
3388 : 8073984 : ffp=zero
3389 : 8073984 : ffpp=zero
3390 [ + + ]: 8073984 : if(ii<=npti-3) then
3391 : 8049168 : ffp=(four*fct(ii+1)-three*fct(ii)-fct(ii+2))/dx/two
3392 : : ffpp=(-fct(ii+3)+four*fct(ii+2)-five*fct(ii+1)+&
3393 : 8049168 : & two*fct(ii))/dx/dx
3394 : : endif
3395 : 8073984 : xx2=xx+half*dx
3396 : 8073984 : xx1=zero
3397 [ + + ]: 8073984 : if(ii>1) xx1=xx-half*dx
3398 : : xsum=xsum+ff*log(abs((xx2-pole)/(xx1-pole)))+ffp*(xx2-xx1+(pole-xx)*log(abs((xx2-pole)/(xx1-pole))))+&
3399 : 8082256 : & half*ffpp*((xx-pole)**2*log(abs((xx2-pole)/(xx1-pole)))+(xx2**2-xx1**2+(two*pole-four*xx)*(xx2-xx1))/two)
3400 : : enddo
3401 [ + + ]: 8272 : if(pole<tol3) xsum=zero
3402 : :
3403 : : ! Calculate the derivated optical quantities and output the value
3404 : 8272 : sigma2=(-two/pi)*xsum
3405 : 8272 : eps1=one-(four_pi*sigma2/(pole))
3406 : 8272 : eps2=four*fcti(ip)*pi/(pole)
3407 : :
3408 : : ! A special treatment of the case where eps2 is very small compared to eps1 is needed
3409 [ + + ]: 8272 : if(eps2**2 > eps1**2 * tol12)then
3410 : 4793 : nomega=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
3411 : 4793 : komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
3412 : 4793 : abso=four_pi*fcti(ip)/(nomega*Speed_Light)
3413 [ + + ]: 3479 : else if(eps1>zero)then
3414 : 3214 : nomega=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
3415 : 3214 : komega=half*abs(eps2/sqrt(eps1))
3416 : 3214 : abso=four_pi*fcti(ip)/(nomega*Speed_Light)
3417 [ + - ]: 265 : else if(eps1<zero)then
3418 : 265 : nomega=half*abs(eps2/sqrt(-eps1))
3419 : 265 : komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
3420 : 265 : abso=two*sqrt(-eps1)*pole/(Speed_Light)
3421 : : end if
3422 : :
3423 : 8272 : epsc=cmplx(eps1,eps2,kind=dp)
3424 : 8272 : cos_phi=cmplx(cos(phi),kind=dp)
3425 : 8272 : sin_phi=cmplx(sin(phi),kind=dp)
3426 : 8272 : sqroot=sqrt(epsc-sin_phi*sin_phi)
3427 : 8272 : crefl_s=(cos_phi-sqroot)/(cos_phi+sqroot)
3428 : 8272 : crefl_p=(epsc*cos_phi-sqroot)/(epsc*cos_phi+sqroot)
3429 : :
3430 : 8272 : refl_s=real(crefl_s*conjg(crefl_s))
3431 : 8272 : refl_p=real(crefl_p*conjg(crefl_p))
3432 : :
3433 : 8272 : sigma1=fcti(ip)
3434 [ + - ]: 8272 : if(au_units==0) then
3435 : 8272 : pole=pole*Ha_eV
3436 : 8272 : sigma1=sigma1*Ohmcm
3437 : 8272 : sigma2=sigma2*Ohmcm
3438 : 8272 : abso=abso*Ohmcm*Speed_Light/(Speed_Light_SI*four_pi*eps0)
3439 : : endif
3440 : 8272 : write(eps_unt,'(5e18.10)') pole,sigma1,sigma2,eps1,eps2
3441 : 8281 : write(abs_unt,'(6e18.10)') pole,nomega,komega,refl_s,refl_p,abso
3442 : :
3443 : : end do
3444 : :
3445 : 9 : close(eps_unt)
3446 : 9 : close(abs_unt)
3447 : :
3448 [ - + ]: 9 : ABI_FREE(fct)
3449 : :
3450 : 9 : end subroutine msig
3451 : : !!***
3452 : :
3453 : : end module m_conducti
3454 : : !!***
|
