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 : module m_conducti
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_xmpi
30 : use m_wfk
31 : use m_hdr
32 : use m_nctk
33 : use netcdf
34 :
35 : use defs_abitypes, only : MPI_type
36 : use m_io_tools, only : open_file, close_unit, get_unit
37 : use m_fstrings, only : sjoin
38 : use m_matrix, only : matr3inv
39 : use m_hide_lapack, only : jacobi
40 : use m_occ, only : getnel
41 : use m_geometry, only : metric
42 : use m_splines, only : splint,spline
43 : use m_mpinfo, only : distrb2,init_mpi_enreg,destroy_mpi_enreg,proc_distrb_cycle,initmpi_band
44 :
45 : implicit none
46 :
47 : private
48 : !!***
49 :
50 : public :: conducti_paw
51 : public :: conducti_paw_core
52 : public :: conducti_nc
53 :
54 : !I/O parameters
55 : !Set to true to use netcdf-MPIIO when available
56 : logical,parameter :: use_netcdf_mpiio=.true.
57 :
58 : !!***
59 :
60 : contains
61 : !!***
62 :
63 : !!****f* m_conducti/conducti_paw
64 : !! NAME
65 : !! conducti_paw
66 : !!
67 : !! FUNCTION
68 : !! This program computes the elements of the optical frequency dependent
69 : !! conductivity tensor and the conductivity along the three principal axes
70 : !! from the Kubo-Greenwood formula for PAW formalism
71 : !!
72 : !! INPUTS
73 : !! filnam=generic name for input data
74 : !! filnam_out=generic name for output data
75 : !! [varocc]=if true, read arbitrary occupations from a file
76 : !!
77 : !! OUTPUT
78 : !! Only printing
79 : !!
80 : !! NOTES
81 : !! bantot
82 : !! doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy.
83 : !! dom=frequency range
84 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
85 : !! eigen11(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
86 : !! in direction x
87 : !! eigen12(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
88 : !! in direction y
89 : !! eigen13(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree)
90 : !! in direction z
91 : !! ecut=kinetic energy planewave cutoff (hartree).
92 : !! fermie= fermi energy (Hartree)
93 : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$).
94 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
95 : !! kin11= Onsager kinetic coeficient=optical conductivity
96 : !! kin12= Onsager kinetic coeficient
97 : !! kin21= Onsager kinetic coeficient
98 : !! kin22= Onsager kinetic coeficient
99 : !! Kth=thermal conductivity
100 : !! mom=number of frequency for conductivity computation
101 : !! mband=maximum number of bands.
102 : !! natom = number of atoms in the unit cell.
103 : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
104 : !! nkpt=number of k points in the IBZ for this perturbation
105 : !! ngfft(3)=integer fft box dimensions.
106 : !! nspinor=number of spinorial components of the wavefunctions.
107 : !! nsppol=1 for unpolarized, 2 for spin-polarized.
108 : !! ntypat = number of atom types.
109 : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
110 : !! occopt==option for occupancies
111 : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor))
112 : !! rprimd(3,3)=real space primitive translations.
113 : !! of primitive translations.
114 : !! Sth=thermopower
115 : !! tsmear=smearing width (or temperature) in Hartree
116 : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
117 : !! wind=frequency windows for computations of sigma
118 : !! wtk(nkpt)=weight assigned to each k point.
119 : !! znucl(natom)=atomic number of atoms
120 : !! np_sum=noziere-pines sumrule
121 : !!
122 : !! SOURCE
123 :
124 2 : subroutine conducti_paw(filnam,filnam_out,varocc)
125 :
126 : !Arguments -----------------------------------
127 : !scalars
128 : character(len=fnlen),intent(in) :: filnam,filnam_out
129 : integer, optional, intent(in) :: varocc
130 :
131 : !Local variables-------------------------------
132 : !scalars
133 : integer,parameter :: master=0
134 : integer :: bsize,bd_stride,dimid,iomode,bantot,bdtot_index,ncid,varid,nb_per_proc,etiq
135 : integer :: comm,fform1,headform,iband,ijband,ierr,ikpt,master_band,idum,iproc,num_tasks_max
136 : integer :: iom,isppol,jband,l1,l2,mband,me,mpierr,mom
137 : integer :: natom,nband_k,nkpt,nproc,nspinor,nsppol,ntypat,broad_mode
138 : integer :: occopt,iunt,opt_unt,occ_unt,iocc,my_iband,pnp_size,add_drude
139 : integer :: lij_unt,sig_unt,sigd_unt,kth_unt,ocond_unt,occunit,occnpt,au_units
140 : logical :: nc_unlimited,mykpt,myband,iomode_estf_mpiio,read_half_dipoles,omega0
141 : real(dp) :: dirac,del,deltae,deltae_min,deltae_min_tmp,dhdk2_g,diff_eig,diff_occ
142 : real(dp) :: dosdeltae,ecut,entropy,fermie,fermih,kin_fact,maxocc,docc_deig
143 : real(dp) :: np_sum,np_sum_k1,np_sum_2,omin,omax,dom,oml,sig,socc,socc_k,fact_omega0
144 : real(dp) :: Tatm,tphysel,tsmear,ucvol,eig_in_max,eig_in_min,phi
145 : character(len=fnlen) :: filnam1,filnam_gen,occfile
146 : character(len=500) :: msg
147 : character(len=100) :: line
148 2 : type(hdr_type) :: hdr
149 2 : type(MPI_type) :: mpi_enreg
150 : !arrays
151 : 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)
152 : integer :: dummy(0,0,0)
153 2 : integer,allocatable :: nband(:),num_tasks(:)
154 : real(dp) :: dhdk2_r(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
155 2 : real(dp),allocatable :: cond_nd(:,:,:),cond_nd_k(:,:,:)
156 2 : real(dp),allocatable :: doccde(:),doccde_k(:),eig0_k(:),eigen0(:),eig0nc(:,:,:)
157 2 : real(dp),allocatable :: occ(:),occ_k(:),wtk(:),oml1(:),occ_in(:),eig_in(:),occ_tmp(:),ypp(:)
158 2 : real(dp),allocatable :: kin11(:,:),kin12(:),kin21(:),kin22(:)
159 2 : real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
160 2 : real(dp),allocatable :: psinablapsi(:,:,:),sig_abs(:)
161 :
162 : ! *********************************************************************************
163 :
164 : ! ---------------------------------------------------------------------------------
165 : ! Read input data
166 :
167 : !Global MPI communicator
168 2 : comm = xmpi_world
169 2 : nproc = xmpi_comm_size(comm)
170 2 : me = xmpi_comm_rank(comm)
171 :
172 : !Read input parameters file
173 2 : if (me==master) then
174 2 : if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then
175 0 : ABI_ERROR(msg)
176 : end if
177 2 : rewind(iunt)
178 2 : read(iunt,*)
179 2 : read(iunt,'(a)') filnam_gen ! Generic name for the files
180 2 : filnam1=trim(filnam_gen)//'_OPT'
181 : ! Read frequency range
182 2 : read(iunt,*) dom,omin,omax,mom
183 2 : read(iunt,end=13,fmt='(a)') line
184 2 : read(line,end=11,fmt=*) broad_mode,au_units,phi,add_drude
185 0 : goto 14
186 2 : 11 read(line,end=13,fmt=*) broad_mode,au_units
187 0 : phi=zero; add_drude=0
188 0 : goto 14
189 2 : 13 broad_mode=1 ; au_units=0; phi=zero; add_drude=0
190 : 14 continue
191 : ! In case of varocc read filename of occupation datafile
192 2 : if (present(varocc)) then
193 0 : if (me==master) then
194 0 : write(std_out,*) 'Warning, this undocumented feature is highly experimental'
195 0 : write(std_out,*) 'and of limited physical validity, proceed with extreme caution!!!'
196 : end if
197 0 : read(iunt,*) occfile !filename of the non-eq distribution function
198 : end if
199 2 : close(iunt)
200 : end if
201 2 : phi=phi*pi/180.0_dp
202 :
203 : !Send data to all procs
204 2 : call xmpi_bcast(dom,master,comm,mpierr)
205 2 : call xmpi_bcast(omin,master,comm,mpierr)
206 2 : call xmpi_bcast(omax,master,comm,mpierr)
207 2 : call xmpi_bcast(mom,master,comm,mpierr)
208 2 : call xmpi_bcast(broad_mode,master,comm,mpierr)
209 2 : call xmpi_bcast(au_units,master,comm,mpierr)
210 2 : call xmpi_bcast(phi,master,comm,mpierr)
211 2 : call xmpi_bcast(add_drude,master,comm,mpierr)
212 :
213 :
214 : ! ---------------------------------------------------------------------------------
215 : ! Read OPT file
216 :
217 : !Check for FORTRAN/.nc OPT file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF
218 2 : if (me==master) then
219 2 : call nctk_fort_or_ncfile(filnam1,iomode,msg)
220 2 : if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
221 : end if
222 2 : call xmpi_bcast(filnam1,master,comm,mpierr)
223 2 : call xmpi_bcast(iomode,master,comm,mpierr)
224 :
225 : !Open OPT file and read HEADER
226 2 : if (me==master) then
227 2 : if (iomode==IO_MODE_ETSF) then
228 2 : NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self))
229 2 : call hdr%ncread(ncid,fform1)
230 : else
231 0 : if (open_file(filnam1,msg,newunit=opt_unt,form="unformatted",status="old")/=0) then
232 0 : ABI_ERROR(msg)
233 : end if
234 0 : call hdr%fort_read(opt_unt,fform1,rewind=.true.)
235 : end if
236 2 : ABI_CHECK(fform1/=0,sjoin("Error while reading ",filnam1))
237 2 : ABI_CHECK(fform1==610.or.fform1==620,"Conducti requires an OPT file with fform=610 or 620!")
238 : end if
239 2 : call hdr%bcast(master,me,comm)
240 2 : call xmpi_bcast(fform1,master,comm,mpierr)
241 :
242 : !Extract info from the header
243 2 : headform=hdr%headform
244 2 : bantot=hdr%bantot
245 2 : ecut=hdr%ecut_eff
246 2 : natom=hdr%natom
247 2 : nkpt=hdr%nkpt
248 2 : nspinor=hdr%nspinor
249 2 : nsppol=hdr%nsppol
250 2 : ntypat=hdr%ntypat
251 2 : occopt=hdr%occopt
252 26 : rprimd(:,:)=hdr%rprimd(:,:)
253 2 : fermie=hdr%fermie
254 2 : tsmear=hdr%tsmear
255 6 : ABI_MALLOC(nband,(nkpt*nsppol))
256 6 : ABI_MALLOC(occ,(bantot))
257 6 : ABI_MALLOC(wtk,(nkpt))
258 4 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
259 32 : occ(1:bantot)=hdr%occ(1:bantot)
260 4 : wtk(1:nkpt)=hdr%wtk(1:nkpt)
261 4 : mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt)
262 2 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell
263 :
264 : !Read eigenvalues
265 6 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
266 2 : if (me==master) then
267 2 : if (iomode==IO_MODE_ETSF) then
268 2 : varid=nctk_idname(ncid,"eigenvalues")
269 10 : ABI_MALLOC(eig0nc,(mband,nkpt,nsppol))
270 2 : NCF_CHECK(nf90_get_var(ncid,varid,eig0nc))
271 4 : eigen0 = reshape(eig0nc,[mband*nkpt*nsppol])
272 2 : ABI_FREE(eig0nc)
273 : !Close file here because the rest will possibly be read with collective I/O
274 2 : NCF_CHECK(nf90_close(ncid))
275 : else
276 0 : read(opt_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
277 : end if
278 : end if
279 2 : call xmpi_bcast(eigen0,master,comm,mpierr)
280 :
281 : !In case of varocc, use arbitrary occupations
282 : !Read occfile and overwrite occ by interpolating the data in the input file
283 2 : if (present(varocc)) then
284 0 : if (me==master) then
285 0 : ABI_MALLOC(occ_tmp,(bantot))
286 0 : if (open_file(occfile,msg,newunit=occ_unt,form='formatted')/=0) then
287 0 : ABI_ERROR(msg)
288 : end if
289 : ! Get units used in occfile (1=Ha, 2= eV)
290 0 : read(occ_unt,*) occunit
291 0 : read(occ_unt,*) occnpt
292 0 : ABI_MALLOC(occ_in,(occnpt))
293 0 : ABI_MALLOC(eig_in,(occnpt))
294 0 : do iocc=1,occnpt
295 0 : read(occ_unt,*) eig_in(iocc),occ_in(iocc)
296 0 : if (occunit==2) then !Convert units from eV to Hartree
297 0 : eig_in(iocc)=eig_in(iocc)/Ha_eV
298 : occ_in(iocc)=occ_in(iocc)
299 : end if
300 : end do
301 0 : close(occ_unt)
302 : ! Interpolation
303 0 : ABI_MALLOC(ypp,(occnpt))
304 0 : call spline(eig_in,occ_in,occnpt,zero,zero,ypp)
305 : ! Interpolate neccessary values
306 0 : eig_in_max=maxval(eig_in)
307 0 : eig_in_min=minval(eig_in)
308 0 : do iocc=1,bantot
309 : ! Check for Extrapolation and set them to physically sound values
310 0 : if (eigen0(iocc)<eig_in_min) then
311 0 : occ_tmp(iocc)=one
312 0 : else if (eigen0(iocc)>eig_in_max) then
313 0 : occ_tmp(iocc)=zero
314 : else
315 0 : call splint(occnpt,eig_in,occ_in,ypp,1,eigen0(iocc),occ_tmp(iocc),ierr)
316 : end if
317 : end do
318 0 : occ=occ_tmp
319 : ! Clean up
320 0 : ABI_FREE(ypp)
321 0 : ABI_FREE(occ_tmp)
322 0 : ABI_FREE(occ_in)
323 0 : ABI_FREE(eig_in)
324 : end if
325 0 : call xmpi_bcast(occ,master,comm,mpierr)
326 0 : occopt=2
327 : end if ! varocc?
328 :
329 : !---------------------------------------------------------------------------------
330 : ! Prepare kpt/band parallelization
331 :
332 2 : call init_mpi_enreg(mpi_enreg)
333 2 : mpi_enreg%comm_kpt=comm
334 2 : mpi_enreg%me_kpt=me
335 2 : mpi_enreg%nproc_spkpt=nproc
336 2 : mpi_enreg%paralbd=1
337 10 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol))
338 6 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt))
339 2 : call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg)
340 2 : call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol)
341 :
342 : !---------------------------------------------------------------------------------
343 : !Print some data
344 :
345 2 : Tatm=tsmear*Ha_K
346 2 : if (me==master) then
347 2 : write(std_out,*)
348 2 : write(std_out,'(a)' )' Input data:'
349 2 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
350 2 : write(std_out,*)
351 2 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
352 2 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
353 2 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
354 2 : write(std_out,'(a,i8)') ' natom =',natom
355 2 : write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol
356 2 : write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
357 2 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
358 2 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
359 : end if
360 :
361 : ! ---------------------------------------------------------------------------------
362 : ! Compute derivative of occupations wrt the energy
363 :
364 6 : ABI_MALLOC(doccde,(mband*nkpt*nsppol))
365 2 : if (occopt<=2) then
366 0 : if (me==master) then
367 0 : write(std_out,'(a,i4)') ' occopt =',occopt
368 : end if
369 0 : doccde=zero
370 : else
371 2 : tphysel=zero
372 2 : maxocc=two/(nsppol*nspinor)
373 : call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermih,maxocc,mband,nband,&
374 2 : & socc,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,12,wtk)
375 2 : entropy=tsmear*entropy
376 2 : if (me==master) then
377 2 : write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha) =',entropy
378 : end if
379 2 : entropy=entropy/socc
380 2 : if (me==master) then
381 2 : write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha/e) =',entropy
382 : end if
383 : endif
384 :
385 : !---------------------------------------------------------------------------------
386 : ! Determine the frequency range and allocate frequency-dependent arrays
387 :
388 2 : del=(omax-omin)/(mom-1)
389 6 : ABI_MALLOC(oml1,(mom))
390 2002 : do iom=1,mom
391 2002 : oml1(iom)=omin+dble(iom-1)*del
392 : end do
393 :
394 8 : ABI_MALLOC(kin11,(mom,nsppol))
395 4 : ABI_MALLOC(kin12,(mom))
396 4 : ABI_MALLOC(kin21,(mom))
397 4 : ABI_MALLOC(kin22,(mom))
398 6 : ABI_MALLOC(cond_nd,(3,3,mom))
399 4 : ABI_MALLOC(sig_abs,(mom))
400 4 : ABI_MALLOC(Kth,(mom))
401 4 : ABI_MALLOC(Stp,(mom))
402 2004 : kin11 = zero
403 2002 : kin12 = zero
404 2002 : kin21 = zero
405 2002 : kin22 = zero
406 26002 : cond_nd = zero
407 2002 : sig_abs = zero
408 2002 : Kth = zero
409 2002 : Stp = zero
410 :
411 : !---------------------------------------------------------------------------------
412 : !Prepare valence-valence dipoles reading
413 :
414 2 : iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
415 :
416 : !In case of netCDF access to OPT file, prepare collective I/O
417 2 : if (iomode == IO_MODE_ETSF) then
418 2 : if (iomode_estf_mpiio) then
419 2 : NCF_CHECK(nctk_open_read(ncid,filnam1,comm))
420 2 : varid=nctk_idname(ncid,"dipole_valence_valence")
421 2 : if (nproc>1) then
422 0 : NCF_CHECK(nctk_set_collective(ncid,varid))
423 0 : ABI_MALLOC(num_tasks,(nproc))
424 0 : do iproc=1,nproc
425 0 : num_tasks(iproc)=count(mpi_enreg%proc_distrb==iproc-1)
426 : enddo
427 0 : num_tasks_max=maxval(num_tasks)
428 : end if
429 2 : nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR)
430 2 : read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR)
431 : else
432 0 : if (me==master) then
433 0 : NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self))
434 0 : varid=nctk_idname(ncid,"dipole_valence_valence")
435 : !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
436 : ! NCF_CHECK(nctk_set_collective(ncid,varid))
437 : !end if
438 0 : nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR)
439 0 : read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR)
440 : end if
441 0 : call xmpi_bcast(nc_unlimited,master,comm,ierr)
442 0 : call xmpi_bcast(read_half_dipoles,master,comm,ierr)
443 : end if
444 2 : if (nc_unlimited.and.read_half_dipoles) then
445 0 : msg="The OPT file has a wrong format!"
446 0 : ABI_BUG(msg)
447 : end if
448 : else
449 0 : read_half_dipoles=(fform1==620)
450 : end if
451 :
452 2 : if (iomode_estf_mpiio) then
453 : !If MPI-IO, store only ib elements for each jb
454 6 : ABI_MALLOC(psinablapsi,(2,3,mband))
455 : else
456 : !If not, store all pairs (or half)
457 0 : if (read_half_dipoles) then ! only ib>=jb
458 0 : ABI_MALLOC(psinablapsi,(2,3,(mband*(mband+1))/2))
459 : else
460 0 : ABI_MALLOC(psinablapsi,(2,3,mband*mband))
461 : end if
462 : end if
463 8 : pnp_size=size(psinablapsi)
464 :
465 : !---------------------------------------------------------------------------------
466 : ! Compute conductivity
467 :
468 6 : ABI_MALLOC(kin11_k,(mom))
469 4 : ABI_MALLOC(kin12_k,(mom))
470 4 : ABI_MALLOC(kin21_k,(mom))
471 4 : ABI_MALLOC(kin22_k,(mom))
472 6 : ABI_MALLOC(cond_nd_k,(3,3,mom))
473 :
474 2 : np_sum = zero
475 2 : socc = zero
476 2 : deltae = zero
477 2 : deltae_min = 1.d99
478 :
479 : !LOOP OVER SPINS/K
480 2 : bdtot_index = 0
481 4 : do isppol=1,nsppol
482 6 : do ikpt=1,nkpt
483 2 : etiq=ikpt+(isppol-1)*nkpt
484 2 : nband_k=nband(ikpt+(isppol-1)*nkpt)
485 2 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))
486 32 : master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
487 :
488 : ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point
489 : ! Master node reads and send to relevant processor
490 2 : if (.not.iomode_estf_mpiio.and.me==master) then
491 0 : if (iomode==IO_MODE_ETSF) then
492 0 : psinablapsi=zero
493 0 : if (nc_unlimited) then
494 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]
495 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
496 0 : else if (.not.read_half_dipoles) then
497 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]
498 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
499 : else
500 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]
501 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
502 : end if
503 : else
504 0 : psinablapsi=zero
505 0 : bsize=nband_k**2;if (read_half_dipoles) bsize=(nband_k*(nband_k+1))/2
506 0 : read(opt_unt)(psinablapsi(1:2,1,ijband),ijband=1,bsize)
507 0 : read(opt_unt)(psinablapsi(1:2,2,ijband),ijband=1,bsize)
508 0 : read(opt_unt)(psinablapsi(1:2,3,ijband),ijband=1,bsize)
509 : end if
510 0 : if (.not.mykpt) then
511 0 : call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,master_band,comm,etiq,ierr)
512 : end if
513 : end if
514 :
515 : ! Select k-points for current proc
516 2 : if (mykpt) then
517 :
518 6 : ABI_MALLOC(eig0_k,(nband_k))
519 4 : ABI_MALLOC(occ_k,(nband_k))
520 4 : ABI_MALLOC(doccde_k,(nband_k))
521 :
522 26002 : cond_nd_k = zero
523 2002 : kin11_k = zero
524 2002 : kin12_k = zero
525 2002 : kin21_k = zero
526 2002 : kin22_k = zero
527 2 : np_sum_k1 = zero
528 2 : socc_k = zero
529 :
530 : ! k-dependent data
531 32 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
532 32 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
533 32 : doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
534 :
535 : ! In case of non MPI-IO, receive all (n,m) dipoles from master proc
536 : ! Then broadcast them to all band processors
537 2 : if (.not.iomode_estf_mpiio) then
538 0 : if (me/=master.and.me==master_band) then
539 0 : call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,me,comm,etiq,ierr)
540 : end if
541 0 : call xmpi_bcast(psinablapsi,master,mpi_enreg%comm_band,mpierr)
542 : end if
543 :
544 : ! LOOP OVER BANDS n
545 32 : do iband=1,nband_k
546 :
547 : !If MPI-IO, store only ib elements for each jb
548 : !If not, store all (ib,jb) pairs
549 30 : my_iband=merge(1,iband,iomode_estf_mpiio)
550 :
551 : ! Select bands for current proc
552 30 : myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me)
553 32 : if (myband) then
554 :
555 : ! In case of MPI-IO, read valence-valence dipoles for band n
556 30 : if (iomode_estf_mpiio) then
557 30 : if (nc_unlimited) then
558 0 : nc_start_6=[1,1,iband,ikpt,isppol,1] ; nc_count_6=[2,3,1,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1]
559 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
560 30 : else if (.not.read_half_dipoles) then
561 0 : nc_start_6=[1,1,1,iband,ikpt,isppol] ; nc_count_6=[2,3,mband,1,1,1] ; nc_stride_6=[1,1,1,1,1,1]
562 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
563 : else
564 330 : nc_start_5=[1,1,(iband*(iband-1))/2+1,ikpt,isppol] ; nc_count_5=[2,3,iband,1,1] ; nc_stride_5=[1,1,1,1,1]
565 30 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
566 : end if
567 : end if
568 :
569 : ! LOOP OVER BANDS m
570 295 : do jband=1,iband
571 265 : diff_occ = occ_k(iband)-occ_k(jband)
572 265 : diff_eig = eig0_k(iband)-eig0_k(jband)
573 :
574 265 : dhdk2_r = zero
575 265 : dhdk2_g = zero
576 :
577 265 : if (read_half_dipoles) then
578 265 : ijband=(my_iband*(my_iband-1))/2+jband
579 : else
580 : !psinablapsi size is mband for netCDF I/O, nband_k for Fortran I/O
581 0 : bd_stride=merge(mband,nband_k,iomode==IO_MODE_ETSF)
582 0 : ijband=(my_iband-1)*bd_stride+jband
583 : end if
584 :
585 1060 : do l2=1,3
586 3445 : do l1=1,3
587 : dhdk2_r(l1,l2)=dhdk2_r(l1,l2)+(&
588 : & psinablapsi(1,l1,ijband)*psinablapsi(1,l2,ijband)&
589 3180 : & +psinablapsi(2,l1,ijband)*psinablapsi(2,l2,ijband))
590 : end do
591 : end do
592 1060 : do l1=1,3
593 : dhdk2_g=dhdk2_g &
594 : & +(psinablapsi(1,l1,ijband)*psinablapsi(1,l1,ijband) &
595 1060 : & +psinablapsi(2,l1,ijband)*psinablapsi(2,l1,ijband))*third ! Average over directions
596 : end do
597 :
598 :
599 : !Minimal validity limit
600 265 : deltae_min_tmp=dabs(diff_eig)
601 265 : if ((deltae_min_tmp>=tol5).and.(deltae_min_tmp<=deltae_min)) deltae_min=deltae_min_tmp
602 :
603 : !Conductivity for each omega - Apply KG formula
604 265 : kin_fact=(eig0_k(iband)+eig0_k(jband))*half-(fermie+entropy)
605 265 : omega0=dabs(diff_occ)<tol12
606 265 : fact_omega0=merge(two,one,omega0)
607 : !Evaluate sumrule
608 : if (.not.omega0) then
609 213 : docc_deig=dabs(diff_occ/diff_eig)
610 52 : else if(add_drude==1) then
611 0 : docc_deig=dabs(doccde_k(iband))
612 : else
613 : docc_deig=zero
614 : endif
615 265 : np_sum_k1=np_sum_k1 + two*dhdk2_g*docc_deig/fact_omega0
616 : !
617 265295 : do iom=1,mom
618 265000 : oml=oml1(iom)
619 265000 : if (.not.omega0) then
620 213000 : docc_deig=dabs(diff_occ/oml)
621 52000 : else if(add_drude==1) then
622 0 : docc_deig=dabs(doccde_k(iband))
623 : else
624 : docc_deig=zero
625 : endif
626 265000 : if(broad_mode==0) then
627 : dirac=dexp(-((abs(diff_eig)-oml)/(sqrt(two)*dom))**2)/(dom*dsqrt(pi*two))-& ! Take into account (n,m) and (m,n)
628 0 : & dexp(-((abs(diff_eig)+oml)/(sqrt(two)*dom))**2)/(dom*dsqrt(pi*two))
629 : else
630 265000 : dirac=dom/((abs(diff_eig)-oml)**2+dom**2)/pi-dom/((abs(diff_eig)+oml)**2+dom**2)/pi
631 : endif
632 265000 : sig=dhdk2_g*docc_deig*dirac*pi/(ucvol)
633 265000 : kin11_k(iom)=kin11_k(iom)+sig
634 265000 : kin12_k(iom)=kin12_k(iom)-sig*kin_fact
635 265000 : kin21_k(iom)=kin21_k(iom)-sig*kin_fact
636 265000 : kin22_k(iom)=kin22_k(iom)+sig*kin_fact**2
637 1060265 : do l2=1,3
638 3445000 : do l1=1,3
639 3180000 : cond_nd_k(l1,l2,iom)=cond_nd_k(l1,l2,iom)+dhdk2_r(l1,l2)*docc_deig*dirac*pi/ucvol
640 : end do
641 : end do
642 : end do
643 :
644 : end do !jband
645 :
646 30 : socc_k=socc_k+occ_k(iband)
647 :
648 : end if ! my band?
649 : end do ! iband
650 :
651 : ! Accumulate k-point contribution
652 2002 : do iom=1,mom
653 2000 : kin11(iom,isppol)=kin11(iom,isppol)+wtk(ikpt)*kin11_k(iom)
654 2000 : kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
655 2000 : kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
656 2000 : kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
657 26002 : cond_nd(:,:,iom)=cond_nd(:,:,iom)+wtk(ikpt)*cond_nd_k(:,:,iom)
658 : end do
659 2 : np_sum=np_sum+wtk(ikpt)*np_sum_k1
660 2 : socc=socc+wtk(ikpt)*socc_k
661 :
662 : ! Validity limit
663 2 : deltae=deltae+(eig0_k(nband_k)-fermie)
664 :
665 2 : ABI_FREE(eig0_k)
666 2 : ABI_FREE(occ_k)
667 2 : ABI_FREE(doccde_k)
668 :
669 : ! End loop over kpt/spin
670 : end if ! My kpt?
671 4 : bdtot_index=bdtot_index+nband_k
672 : end do ! ikpt
673 : end do ! isppol
674 :
675 2 : if (iomode == IO_MODE_ETSF) then
676 2 : if(iomode_estf_mpiio.and.nproc>1) then
677 0 : do idum=num_tasks(me+1)+1,num_tasks_max
678 0 : if (nc_unlimited) then
679 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]
680 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
681 0 : else if (.not.read_half_dipoles) then
682 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]
683 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
684 : else
685 0 : nc_start_5=[1,1,1,1,1] ; nc_count_5=[0,0,0,0,0] ; nc_stride_5=[1,1,1,1,1]
686 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
687 : end if
688 : enddo
689 0 : ABI_FREE(num_tasks)
690 : endif
691 : endif
692 :
693 : !Accumulate kpt/band contributions over processors
694 2 : call xmpi_sum(kin11,comm,mpierr)
695 2 : call xmpi_sum(kin12,comm,mpierr)
696 2 : call xmpi_sum(kin21,comm,mpierr)
697 2 : call xmpi_sum(kin22,comm,mpierr)
698 2 : call xmpi_sum(np_sum,comm,mpierr)
699 2 : call xmpi_sum(socc,comm,mpierr)
700 2 : call xmpi_sum(deltae,comm,mpierr)
701 2 : call xmpi_min(deltae_min,comm,mpierr)
702 2 : call xmpi_sum(cond_nd,comm,mpierr)
703 2 : deltae=deltae/mpi_enreg%nproc_band
704 :
705 2 : ABI_FREE(psinablapsi)
706 :
707 : !---------------------------------------------------------------------------------
708 : ! Output results
709 :
710 : !Print file headers (only master node)
711 2 : if (me==master) then
712 : ! Standard output
713 : write(std_out,'(a,f10.5,a,f10.5,a)')&
714 2 : & ' Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha', &
715 4 : & deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
716 : write(std_out,'(a,f10.5,a,f10.5,a)')&
717 2 : & ' DeltaE min =',deltae_min,' Ha',deltae_min*Ha_eV,' eV'
718 :
719 : ! _Lij file
720 2 : if (open_file(trim(filnam_out)//'_Lij',msg, newunit=lij_unt, form='formatted', action="write") /= 0) then
721 0 : ABI_ERROR(msg)
722 : end if
723 2 : write(lij_unt,'(a)')' # omega(ua) L11 L12 L22'
724 :
725 : ! _sig_up_dn file
726 2 : if(nsppol==2) then
727 0 : if (open_file(trim(filnam_out)//'_sig_up_dn', msg, newunit=sig_unt, form='formatted', action="write") /= 0) then
728 0 : ABI_ERROR(msg)
729 : end if
730 0 : if(au_units>0) then
731 0 : write(sig_unt,'(2a)')' # omega(ua) cond(ua) ',&
732 0 : & ' cond(ua) UP cond(ua) DN'
733 : else
734 0 : write(sig_unt,'(2a)')' # hbar*omega(eV) cond(ohm.cm)-1',&
735 0 : & ' cond(ohm.cm)-1 UP cond(ohm.cm)-1 DN'
736 : endif
737 : end if
738 :
739 : ! sig_dir file
740 2 : if (open_file(trim(filnam_out)//'_sig_tensor', msg, newunit=sigd_unt,form='formatted', action="write") /= 0) then
741 0 : ABI_ERROR(msg)
742 : end if
743 2 : if(au_units>0) then
744 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)'
745 : else
746 : write(sigd_unt,'(a)')' # hbar*omega(eV) cond_xx(ohm.cm)-1 cond_yy(ohm.cm)-1 cond_zz(ohm.cm)-1 &
747 2 : & cond_xy(ohm.cm)-1 cond_xz(ohm.cm)-1 cond_yz(ohm.cm)-1'
748 : endif
749 :
750 : ! _Kth file
751 2 : if (open_file(trim(filnam_out)//'_Kth', msg, newunit=kth_unt, form='formatted', action="write") /=0) then
752 0 : ABI_ERROR(msg)
753 : end if
754 : write(kth_unt,'(a)')&
755 2 : & " #Thermal conductivity following B. Holst et al Phys. Rev. B 83 (2011) 235120"
756 2 : if(au_units>0) then
757 0 : write(kth_unt,'(a)')' # omega(ua) thermal cond(ua) thermopower(ua)'
758 : else
759 2 : write(kth_unt,'(a)')' # hbar*omega(eV) thermal cond (W/m/K) thermopower(microohm/K)'
760 : endif
761 :
762 : ! Output file
763 2 : if (open_file(trim(filnam_out)//'.out', msg, newunit=ocond_unt, form='formatted', action="write") /= 0) then
764 0 : ABI_ERROR(msg)
765 : end if
766 2 : write(ocond_unt,'(a)' )'#Conducti output file:'
767 2 : write(ocond_unt,'(a)' )'#Contains all results produced by conducti utility'
768 2 : write(ocond_unt,'(a)' )'# '
769 2 : write(ocond_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =' ,mom,omin,omax,dom,' Ha'
770 2 : write(ocond_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1)
771 2 : write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2)
772 2 : write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3)
773 2 : write(ocond_unt,'(a,i8)' ) '# natom =',natom
774 2 : write(ocond_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol
775 2 : write(ocond_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha'
776 2 : write(ocond_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV'
777 :
778 2 : write(ocond_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin'
779 2 : write(ocond_unt,'(a,f15.5)' )'# Number of electrons = ',socc
780 2 : write(ocond_unt,'(a,f15.5)' )'# sumrule = ',np_sum/socc
781 : write(ocond_unt,'(a,f10.5,a,f10.5,a)' )&
782 2 : & '# Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
783 2 : write(ocond_unt,'(a)' )'# '
784 2 : if(au_units>0) then
785 0 : write(ocond_unt,'(a)')'# omega(ua) cond(ua) thermal cond(ua) thermopower(ua)'
786 : else
787 2 : write(ocond_unt,'(a)')'# hbar*omega(eV) cond(ohm.cm)-1 thermal cond(W/m/K) thermopower(microvolt/K) '
788 : endif
789 :
790 : end if ! me==master?
791 :
792 : !Compute (and print) thermal conductivity and thermopower
793 2002 : do iom=1,mom
794 2000 : oml=oml1(iom)
795 :
796 4000 : do isppol=1,nsppol
797 2000 : if (dabs(kin11(iom,isppol))<tol19) kin11(iom,isppol)=zero
798 4000 : sig_abs(iom)=sig_abs(iom)+kin11(iom,isppol)
799 : end do
800 :
801 2000 : Kth(iom)=kin22(iom)/tsmear
802 2000 : Stp(iom)=zero
803 2000 : if (sig_abs(iom)/=zero) then
804 2000 : Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/sig_abs(iom))/tsmear
805 2000 : Stp(iom)=kin12(iom)/(sig_abs(iom)*tsmear)
806 : else
807 0 : Kth(iom)=Kth(iom-1)
808 0 : Stp(iom)=Stp(iom-1)
809 : endif
810 2000 : if (dabs(Kth(iom))<tol19) Kth(iom)=zero
811 2000 : if (dabs(Stp(iom))<tol19) Stp(iom)=zero
812 2000 : if (abs(kin12(iom))<10.0_dp**(-80)) kin12(iom)=zero
813 2000 : if (abs(kin21(iom))<10.0_dp**(-80)) kin21(iom)=zero
814 2002 : if (abs(kin22(iom))<10.0_dp**(-80)) kin22(iom)=zero
815 : end do
816 :
817 : !Compute the imaginary part of the conductivity (principal value)
818 : ! +derived optical properties.
819 2 : if (me==master) then
820 2 : call msig(sig_abs,mom,oml1,filnam_out,phi,au_units)
821 : end if
822 :
823 : ! Units
824 2 : if (me==master) then
825 2002 : do iom=1,mom
826 2000 : oml=oml1(iom)
827 2002 : write(lij_unt,'(f12.5,3es22.12)') oml,sig_abs(iom),kin12(iom),kin22(iom)
828 : enddo
829 : endif
830 :
831 2 : np_sum_2=zero
832 2000 : do iom=1,mom-1
833 2000 : np_sum_2=np_sum_2+(sig_abs(iom)+sig_abs(iom+1))*(oml1(iom+1)-oml1(iom))/2.0_dp
834 : enddo
835 2 : np_sum_2=2.0_dp*np_sum_2*ucvol/pi
836 2002 : if(au_units==0) oml1=oml1*Ha_eV
837 26002 : if(au_units==0) cond_nd=cond_nd*Ohmcm
838 2002 : if(au_units==0) sig_abs=sig_abs*Ohmcm
839 2004 : if(au_units==0) kin11=kin11*Ohmcm
840 2002 : if(au_units==0) Kth=Kth*3.4057d9/Ha_K
841 2002 : if(au_units==0) Stp=Stp*Ha_J/(Ha_K*e_CB)*10.0_dp**6
842 :
843 :
844 2 : if (me==master) then
845 2 : write(std_out,'(a,f15.5)' )'# Number of electrons = ',socc
846 2 : write(std_out,'(a,f15.5)' )'# sumrule = ',np_sum/socc
847 2 : write(std_out,'(a,f15.5)' )'# sumrule (integration) = ',np_sum_2/socc
848 2002 : do iom=1,mom
849 2000 : oml=oml1(iom)
850 2000 : write(sigd_unt,'(f12.5,6es22.12)') oml,cond_nd(1,1,iom),cond_nd(2,2,iom),cond_nd(3,3,iom),&
851 4000 : & cond_nd(1,2,iom),cond_nd(1,3,iom),cond_nd(2,3,iom)
852 2000 : if (nsppol==2) then
853 0 : write(sig_unt,'(f12.5,3es22.12)') oml,sig_abs(iom), &
854 0 : & kin11(iom,1),kin11(iom,2)
855 : end if
856 2000 : write(kth_unt,'(f12.5,3es22.12)') oml,Kth(iom),Stp(iom)
857 2002 : write(ocond_unt,'(1f12.5,3es22.12)') oml,sig_abs(iom),Kth(iom),Stp(iom)
858 : enddo
859 : end if
860 :
861 :
862 : !---------------------------------------------------------------------------------
863 : ! End
864 :
865 : !Close all files
866 2 : if (me==master) then
867 2 : write(std_out,'(2a)')ch10,'OUTPUT'
868 2 : write(std_out,'(a)')trim(filnam_out)//'_Lij : Onsager kinetic coefficients'
869 2 : write(std_out,'(a)')trim(filnam_out)//'_eps : Optical conductivity and dielectric function'
870 2 : write(std_out,'(a)')trim(filnam_out)//'_sig_tensor : Optical conductivity tensor'
871 2 : write(std_out,'(a)')trim(filnam_out)//'_Kth : Thermal conductivity and thermopower'
872 2 : write(std_out,'(a)')trim(filnam_out)//'_abs : n, k, reflectivity, absorption'
873 2 : close(lij_unt)
874 2 : if(nsppol==2) close(sig_unt)
875 2 : close(kth_unt)
876 2 : close(ocond_unt)
877 2 : close(sigd_unt)
878 : end if
879 2 : if (iomode == IO_MODE_ETSF) then
880 2 : if (iomode_estf_mpiio.or.me==master) then
881 2 : NCF_CHECK(nf90_close(ncid))
882 : end if
883 0 : else if (me==master) then
884 0 : ierr=close_unit(opt_unt,msg)
885 0 : ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam1))
886 : end if
887 :
888 : !Release memory space
889 2 : ABI_FREE(kin11)
890 2 : ABI_FREE(kin22)
891 2 : ABI_FREE(kin12)
892 2 : ABI_FREE(kin21)
893 2 : ABI_FREE(kin11_k)
894 2 : ABI_FREE(kin22_k)
895 2 : ABI_FREE(kin12_k)
896 2 : ABI_FREE(kin21_k)
897 2 : ABI_FREE(Stp)
898 2 : ABI_FREE(Kth)
899 2 : ABI_FREE(cond_nd)
900 2 : ABI_FREE(cond_nd_k)
901 2 : ABI_FREE(sig_abs)
902 2 : ABI_FREE(eigen0)
903 2 : ABI_FREE(nband)
904 2 : ABI_FREE(oml1)
905 2 : ABI_FREE(occ)
906 2 : ABI_FREE(doccde)
907 2 : ABI_FREE(wtk)
908 2 : call hdr%free()
909 2 : call destroy_mpi_enreg(mpi_enreg)
910 :
911 4 : end subroutine conducti_paw
912 : !!***
913 :
914 : !----------------------------------------------------------------------
915 :
916 : !!****f* m_conducti/conducti_paw_core
917 : !! NAME
918 : !! conducti_paw_core
919 : !!
920 : !! FUNCTION
921 : !! This program computes the elements of the optical frequency dependent
922 : !! conductivity tensor and the conductivity along the three principal axes
923 : !! from the Kubo-Greenwood formula for PAW formalism
924 : !!
925 : !! INPUTS
926 : !! filnam=generic name for input data
927 : !! filnam_out=generic name for output data
928 : !! [with_absorption]=optiona lflag to activate the computation of absorption (_sigX file) (default=TRUE)
929 : !! [with_emissivity]=optiona lflag to activate the computation of emissivity (_emisX file) (default=FALSE)
930 : !!
931 : !! OUTPUT
932 : !! Only printing
933 : !!
934 : !! NOTES
935 : !! bantot
936 : !! dom=frequency range
937 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
938 : !! ecut=kinetic energy planewave cutoff (hartree).
939 : !! fermie= fermi energy (Hartree)
940 : !! mom=number of frequency for conductivity computation
941 : !! mband=maximum number of bands.
942 : !! natom = number of atoms in the unit cell.
943 : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
944 : !! nkpt=number of k points in the IBZ for this perturbation
945 : !! ngfft(3)=integer fft box dimensions.
946 : !! nspinor=number of spinorial components of the wavefunctions.
947 : !! nsppol=1 for unpolarized, 2 for spin-polarized.
948 : !! ntypat = number of atom types.
949 : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
950 : !! occopt==option for occupancies
951 : !! psinablapsi2(2,3,mband,nphicor,natom)Matrix elements = <Phi_core|Nabla|Phi_i>
952 : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor))
953 : !! rprimd(3,3)=real space primitive translations.
954 : !! of primitive translations.
955 : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
956 : !! wind=frequency window for computations of sigma
957 : !! wtk(nkpt)=weight assigned to each k point.
958 : !!
959 : !! SOURCE
960 :
961 3 : subroutine conducti_paw_core(filnam,filnam_out,with_absorption,with_emissivity)
962 :
963 : !Arguments -----------------------------------
964 : !scalars
965 : character(len=fnlen),intent(in) :: filnam,filnam_out
966 : logical,intent(in),optional :: with_absorption,with_emissivity
967 : !Local variables-------------------------------
968 : !scalars
969 : integer,parameter :: master=0
970 : integer :: iomode,atnbr,bantot,bdtot_index,comm,my_iband,input_atm,itask,num_tasks_max,au_units
971 : integer :: fform2,headform,iatom,iband,icor,ierr,ikpt,iatom_atnbr,itypat,itypat_atnbr
972 : integer :: iom,isppol,l1,mband,me,mom,mpierr,j2,etiq,pnp_size,iproc,broad_mode,absx_unt
973 : integer :: natom,nband_k,nkpt,nphicor,nproc,nspinor,nsppol,ntypat,nphicor_max,natom_atnbr
974 : integer :: occopt,iunt,opt2_unt,ncid,varid,master_band,nb_per_proc,idum
975 : integer :: sigx1_unt,sigx1_up_unt,sigx1_dn_unt,ems_unt,ems_up_unt,ems_dn_unt
976 : logical :: iomode_estf_mpiio,myband,mykpt,need_absorption,need_emissivity
977 : real(dp) :: del_sig,del_emis,deltae,diff_occ,ecut,fermie
978 : real(dp) :: omin,omax,omin_sig,omax_sig,omin_emis,omax_emis
979 : real(dp) :: oml,dom,dom_ctr,dom_max,dom_tan1,dom_tan2,docc_deig
980 : real(dp) :: Tatm,tsmear,ucvol,dirac,diff_eig
981 : character(len=fnlen) :: filnam2,filnam_gen
982 : character(len=500) :: msg
983 : character(len=100) :: line
984 : character(len=20) :: str_atm
985 3 : type(hdr_type) :: hdr
986 3 : type(MPI_type) :: mpi_enreg
987 : !arrays
988 : integer :: nc_count(6),nc_start(6),nc_stride(6)
989 3 : integer,allocatable :: nband(:),ncor(:,:),lcor(:,:),kappacor(:,:),typat(:),num_tasks(:)
990 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),dummy(0,0,0,0,0)
991 3 : real(dp),allocatable :: dom_var1(:,:),dhdk2_g(:)
992 3 : real(dp),allocatable :: eig0_k(:),eigen0(:),eig0nc(:,:,:)
993 3 : real(dp),allocatable :: energy_cor(:,:),edge(:),nphicor_arr(:),occ_cor(:,:),maxocc_cor(:,:)
994 3 : real(dp),allocatable :: occ(:),occ_k(:),wtk(:)
995 3 : real(dp),allocatable :: oml_edge(:,:),oml_emis(:,:)
996 3 : real(dp),allocatable :: psinablapsi2(:,:,:,:)
997 3 : real(dp),allocatable :: sigx1(:,:,:,:),sigx1_av(:,:,:),sigx1_k(:,:,:)
998 3 : real(dp),allocatable :: sum_spin_sigx1(:,:,:),sum_spin_sigx1_av(:,:)
999 3 : real(dp),allocatable :: emisx(:,:,:,:),emisx_av(:,:,:),emisx_k(:,:,:)
1000 3 : real(dp),allocatable :: sum_spin_emisx(:,:,:),sum_spin_emisx_av(:,:)
1001 :
1002 : ! *********************************************************************************
1003 : !optional flags
1004 3 : need_absorption=.true. ;if (present(with_absorption)) need_absorption=with_absorption
1005 3 : need_emissivity=.false.;if (present(with_emissivity)) need_emissivity=with_emissivity
1006 3 : if ((.not.need_absorption).and.(.not.need_emissivity)) return
1007 :
1008 : ! ---------------------------------------------------------------------------------
1009 : ! Read input data
1010 :
1011 : !Global MPI communicators
1012 3 : comm = xmpi_world
1013 3 : nproc = xmpi_comm_size(comm)
1014 3 : me = xmpi_comm_rank(comm)
1015 :
1016 : !Read input parameters file
1017 3 : if (me==master) then
1018 3 : if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then
1019 0 : ABI_ERROR(msg)
1020 : end if
1021 3 : rewind(iunt)
1022 3 : read(iunt,*)
1023 3 : read(iunt,'(a)')filnam_gen ! Generic name for the files
1024 3 : filnam2=trim(filnam_gen)//'_OPT2'
1025 : ! Read frequency range
1026 3 : if (need_absorption) then
1027 3 : read(iunt,fmt='(a)') line
1028 3 : read(line,end=11,fmt=*) dom,omin,omax,mom,input_atm,dom_max,dom_ctr
1029 2 : goto 12
1030 1 : 11 read(line,fmt=*) dom,omin,omax,mom,input_atm ; dom_max=zero ; dom_ctr=zero
1031 : 12 continue
1032 0 : else if (need_emissivity) then
1033 0 : read(iunt,*) dom,omin,omax,mom,input_atm
1034 0 : dom_max=zero;dom_ctr=zero
1035 : end if
1036 3 : read(iunt,end=13,err=13,fmt=*) broad_mode,au_units
1037 1 : goto 14
1038 2 : 13 broad_mode=1 ; au_units=0
1039 : 14 continue
1040 3 : close(iunt)
1041 3 : if (abs(dom_max)>tol10.and.dom_max<dom) then
1042 0 : msg = 'dom_max must be higher than dom!'
1043 0 : ABI_ERROR(msg)
1044 : end if
1045 : end if
1046 :
1047 : !Send data to all procs
1048 3 : call xmpi_bcast(dom,master,comm,mpierr)
1049 3 : call xmpi_bcast(omin,master,comm,mpierr)
1050 3 : call xmpi_bcast(omax,master,comm,mpierr)
1051 3 : call xmpi_bcast(mom,master,comm,mpierr)
1052 3 : call xmpi_bcast(input_atm,master,comm,mpierr)
1053 3 : call xmpi_bcast(dom_max,master,comm,mpierr)
1054 3 : call xmpi_bcast(dom_ctr,master,comm,mpierr)
1055 3 : call xmpi_bcast(broad_mode,master,comm,mpierr)
1056 3 : call xmpi_bcast(au_units,master,comm,mpierr)
1057 :
1058 :
1059 : ! ---------------------------------------------------------------------------------
1060 : ! Read OPT2 file
1061 :
1062 : !Check for FORTRAN/.nc OPT2 file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF
1063 3 : if (me==master) then
1064 3 : call nctk_fort_or_ncfile(filnam2,iomode,msg)
1065 3 : if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
1066 : end if
1067 3 : call xmpi_bcast(filnam2,master,comm,mpierr)
1068 3 : call xmpi_bcast(iomode,master,comm,mpierr)
1069 :
1070 : !Open OPT2 file and read HEADER
1071 3 : if (me==master) then
1072 3 : if (iomode==IO_MODE_ETSF) then
1073 3 : NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self))
1074 3 : call hdr%ncread(ncid,fform2)
1075 : else
1076 0 : if (open_file(filnam2,msg,newunit=opt2_unt,form="unformatted",status="old")/=0) then
1077 0 : ABI_ERROR(msg)
1078 : end if
1079 0 : call hdr%fort_read(opt2_unt,fform2,rewind=.true.)
1080 : end if
1081 3 : ABI_CHECK(fform2/=0,sjoin("Error while reading ",filnam2))
1082 3 : ABI_CHECK(fform2==611.or.fform2==612.or.fform2==613,"OPT2 file format should be fform=611/612/613!")
1083 : end if
1084 3 : call hdr%bcast(master,me,comm)
1085 3 : call xmpi_bcast(fform2,master,comm,mpierr)
1086 :
1087 : !Extract info from the header
1088 3 : headform=hdr%headform
1089 3 : bantot=hdr%bantot
1090 3 : ecut=hdr%ecut_eff
1091 3 : natom=hdr%natom
1092 3 : nkpt=hdr%nkpt
1093 3 : nspinor=hdr%nspinor
1094 3 : nsppol=hdr%nsppol
1095 3 : ntypat=hdr%ntypat
1096 9 : ABI_MALLOC(typat,(natom))
1097 15 : typat = hdr%typat
1098 : occopt=hdr%occopt
1099 39 : rprimd(:,:)=hdr%rprimd(:,:)
1100 3 : fermie=hdr%fermie
1101 3 : tsmear=hdr%tsmear
1102 9 : ABI_MALLOC(nband,(nkpt*nsppol))
1103 9 : ABI_MALLOC(occ,(bantot))
1104 9 : ABI_MALLOC(wtk,(nkpt))
1105 53 : occ(1:bantot)=hdr%occ(1:bantot)
1106 6 : wtk(1:nkpt)=hdr%wtk(1:nkpt)
1107 6 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1108 6 : mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt)
1109 3 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell
1110 3 : itypat_atnbr=typat(input_atm)
1111 3 : natom_atnbr=0
1112 12 : do iatom=1,natom
1113 9 : if(typat(iatom)==itypat_atnbr) natom_atnbr=natom_atnbr+1
1114 12 : if(iatom==input_atm) atnbr=natom_atnbr
1115 : enddo
1116 :
1117 : !Read eigenvalues
1118 9 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
1119 3 : if (me==master) then
1120 3 : if (iomode==IO_MODE_ETSF) then
1121 3 : varid=nctk_idname(ncid,"eigenvalues")
1122 15 : ABI_MALLOC(eig0nc,(mband,nkpt,nsppol))
1123 3 : NCF_CHECK(nf90_get_var(ncid,varid,eig0nc))
1124 6 : eigen0 = reshape(eig0nc,[mband*nkpt*nsppol])
1125 3 : ABI_FREE(eig0nc)
1126 : else
1127 0 : read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
1128 : end if
1129 : end if
1130 3 : call xmpi_bcast(eigen0,master,comm,mpierr)
1131 :
1132 : !Read core states
1133 3 : if (me==master) then
1134 3 : if (iomode==IO_MODE_ETSF) then
1135 3 : NCF_CHECK(nctk_get_dim(ncid,"max_number_of_core_states",nphicor_max))
1136 : else
1137 0 : read(opt2_unt) nphicor_max
1138 : end if
1139 : end if
1140 3 : call xmpi_bcast(nphicor_max,master,comm,mpierr)
1141 :
1142 12 : ABI_MALLOC(ncor,(nphicor_max,ntypat))
1143 12 : ABI_MALLOC(occ_cor,(nphicor_max,ntypat))
1144 9 : ABI_MALLOC(maxocc_cor,(nphicor_max,ntypat))
1145 9 : ABI_MALLOC(lcor,(nphicor_max,ntypat))
1146 9 : ABI_MALLOC(kappacor,(nphicor_max,ntypat))
1147 9 : ABI_MALLOC(energy_cor,(nphicor_max,ntypat))
1148 :
1149 3 : if (me==master) then
1150 9 : ABI_MALLOC(nphicor_arr,(ntypat))
1151 3 : if (iomode==IO_MODE_ETSF) then
1152 3 : varid=nctk_idname(ncid,"n_quantum_number_core")
1153 3 : NCF_CHECK(nf90_get_var(ncid,varid,ncor))
1154 3 : varid=nctk_idname(ncid,"l_quantum_number_core")
1155 3 : NCF_CHECK(nf90_get_var(ncid,varid,lcor))
1156 3 : varid=nctk_idname(ncid,"kappa_core")
1157 3 : NCF_CHECK(nf90_get_var(ncid,varid,kappacor))
1158 3 : varid=nctk_idname(ncid,"eigenvalues_core")
1159 3 : NCF_CHECK(nf90_get_var(ncid,varid,energy_cor))
1160 3 : varid=nctk_idname(ncid,"occupation_core")
1161 3 : NCF_CHECK(nf90_get_var(ncid,varid,occ_cor))
1162 3 : varid=nctk_idname(ncid,"max_occupation_core")
1163 3 : NCF_CHECK(nf90_get_var(ncid,varid,maxocc_cor))
1164 3 : varid=nctk_idname(ncid,"number_of_core_states")
1165 3 : NCF_CHECK(nf90_get_var(ncid,varid,nphicor_arr))
1166 : !Close here netcdf file here because the rest has to be read with collective I/O
1167 3 : NCF_CHECK(nf90_close(ncid))
1168 : else
1169 0 : do itypat=1,ntypat
1170 0 : read(unit=opt2_unt) nphicor_arr(itypat)
1171 0 : do icor=1,nphicor_max
1172 0 : read(unit=opt2_unt,fmt='(a)') line
1173 0 : read(line,end=23,fmt=*) ncor(icor,itypat),lcor(icor,itypat),kappacor(icor,itypat),occ_cor(icor,itypat),maxocc_cor(icor,itypat),&
1174 0 : energy_cor(icor,itypat)
1175 0 : goto 24
1176 0 : 23 read(line,fmt=*) ncor(icor,itypat),lcor(icor,itypat),energy_cor(icor,itypat)
1177 0 : kappacor(icor,itypat)=0
1178 0 : occ_cor(icor,itypat)=one
1179 0 : maxocc_cor(icor,itypat)=one
1180 0 : 24 continue
1181 : enddo
1182 : end do
1183 : end if
1184 3 : nphicor=nphicor_arr(itypat_atnbr)
1185 3 : ABI_FREE(nphicor_arr)
1186 : end if ! master
1187 3 : call xmpi_bcast(ncor,master,comm,mpierr)
1188 3 : call xmpi_bcast(lcor,master,comm,mpierr)
1189 3 : call xmpi_bcast(kappacor,master,comm,mpierr)
1190 3 : call xmpi_bcast(energy_cor,master,comm,mpierr)
1191 3 : call xmpi_bcast(occ_cor,master,comm,mpierr)
1192 3 : call xmpi_bcast(maxocc_cor,master,comm,mpierr)
1193 3 : call xmpi_bcast(nphicor,master,comm,mpierr)
1194 :
1195 9 : ABI_MALLOC(edge,(nphicor))
1196 14 : edge(1:nphicor)=fermie-energy_cor(1:nphicor,itypat_atnbr)
1197 : !---------------------------------------------------------------------------------
1198 : ! Prepare kpt/band parallelization
1199 :
1200 3 : call init_mpi_enreg(mpi_enreg)
1201 3 : mpi_enreg%comm_kpt=comm
1202 3 : mpi_enreg%me_kpt=me
1203 3 : mpi_enreg%nproc_spkpt=nproc
1204 3 : mpi_enreg%paralbd=1
1205 15 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol))
1206 9 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt))
1207 3 : call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg)
1208 3 : call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol)
1209 : !---------------------------------------------------------------------------------
1210 : !Print some data
1211 3 : Tatm=tsmear*Ha_K
1212 3 : if (me==master) then
1213 3 : write(std_out,*)
1214 3 : write(std_out,'(a)')'--------------------------------------------'
1215 3 : write(std_out,'(a,i4)') 'selected atom for X ray emission',input_atm
1216 3 : write(std_out,'(a)')'--------------------------------------------'
1217 3 : if (need_absorption) then
1218 3 : if (abs(dom_max)>tol10) then
1219 1 : write(std_out,'(a)')'************************ ARCTAN SMEARING'
1220 1 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1221 1 : write(std_out,'(a,2f10.5,a)')' dom_max,center =',dom_max,dom_ctr,' Ha'
1222 : else
1223 2 : write(std_out,'(a)')'************************ FIXED SMEARING'
1224 2 : write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1225 : endif
1226 : end if
1227 3 : write(std_out,*)
1228 3 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1,1:3)
1229 3 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(2,1:3)
1230 3 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(3,1:3)
1231 3 : write(std_out,'(a,i8)') ' natom =',natom
1232 3 : write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol
1233 3 : write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
1234 3 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
1235 3 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
1236 3 : write(std_out,*)
1237 3 : write(std_out,*)
1238 3 : write(std_out,'(a)')'--------------------------------------------'
1239 3 : write(std_out,'(a,i4)') ' Number of core orbitals nc=',nphicor
1240 14 : do icor=1,nphicor
1241 14 : if (kappacor(icor,itypat_atnbr)==0) then
1242 3 : write(std_out,'(a,2i4,5f15.5)') ' n, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ', &
1243 3 : ncor(icor,itypat_atnbr),lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1244 3 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1245 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1246 : else
1247 8 : if (kappacor(icor,itypat_atnbr)>0) then
1248 2 : j2=2*lcor(icor,itypat_atnbr)-1
1249 2 : write(std_out,'(a,i4,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1250 2 : & ncor(icor,itypat_atnbr),j2,' / 2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1251 2 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1252 4 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1253 : else
1254 6 : if (kappacor(icor,itypat_atnbr)<-1) then
1255 2 : j2=2*lcor(icor,itypat_atnbr)+1
1256 2 : write(std_out,'(a,i4,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1257 2 : & ncor(icor,itypat_atnbr),j2,'/2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1258 2 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1259 4 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1260 : else
1261 4 : write(std_out,'(a,i4,a,i4,5f15.5)') ' n, j, l, occ, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1262 4 : & ncor(icor,itypat_atnbr),' 1/2',lcor(icor,itypat_atnbr),occ_cor(icor,itypat_atnbr),&
1263 4 : & energy_cor(icor,itypat_atnbr),edge(icor),&
1264 8 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1265 : end if
1266 : end if
1267 : end if
1268 : end do
1269 3 : write(std_out,'(a)')'--------------------------------------------'
1270 : end if ! master
1271 :
1272 : !---------------------------------------------------------------------------------
1273 : ! Determine the frequency range and allocate frequency-dependent arrays
1274 3 : if (need_absorption) then
1275 12 : ABI_MALLOC(oml_edge,(nphicor,mom))
1276 9 : ABI_MALLOC(dom_var1,(nphicor,mom))
1277 3 : omax_sig=omax ; omin_sig=omin
1278 3 : del_sig=(omax_sig-omin_sig)/(mom-1)
1279 3003 : do iom=1,mom
1280 14003 : do icor=1,nphicor
1281 11000 : oml_edge(icor,iom)=-energy_cor(icor,itypat_atnbr)+ dble(iom-1)*del_sig + omin_sig -one
1282 14000 : if ((oml_edge(icor,iom)<=edge(icor)).or.(abs(dom_max)<=tol10)) then
1283 8048 : dom_var1(icor,iom)= dom
1284 : else
1285 2952 : dom_tan1= (oml_edge(icor,iom)-edge(icor))/dom_ctr
1286 2952 : dom_tan2=dom_tan1-one/dom_tan1**2
1287 2952 : dom_var1(icor,iom)=dom+dom_max*(half+piinv*datan(dom_tan2))
1288 : endif
1289 : enddo
1290 : enddo
1291 18 : ABI_MALLOC(sigx1,(nphicor,mom,natom_atnbr,nsppol))
1292 41015 : sigx1=zero
1293 : end if
1294 :
1295 3 : if (need_emissivity) then
1296 4 : ABI_MALLOC(oml_emis,(nphicor,mom))
1297 12 : omin_emis=minval(eigen0)
1298 12 : omax_emis=maxval(eigen0)
1299 1 : del_emis=(omax_emis-omin_emis)/(mom-1)
1300 1001 : do iom=1,mom
1301 4001 : do icor=1,nphicor
1302 4000 : oml_emis(icor,iom)=omin_emis+dble(iom-1)*del_emis!-energy_cor(icor,itypat_atnbr)
1303 : enddo
1304 : end do
1305 6 : ABI_MALLOC(emisx,(nphicor,mom,natom_atnbr,nsppol))
1306 16006 : emisx=zero
1307 : end if
1308 :
1309 : !---------------------------------------------------------------------------------
1310 : ! Prepare core-valence dipoles reading
1311 :
1312 3 : iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
1313 :
1314 : !In case of netCDF access to OPT file, prepare collective I/O
1315 3 : if (iomode == IO_MODE_ETSF) then
1316 3 : if (iomode_estf_mpiio) then
1317 3 : NCF_CHECK(nctk_open_read(ncid,filnam2,comm))
1318 3 : varid=nctk_idname(ncid,"dipole_core_valence")
1319 3 : if (nproc>1) then
1320 0 : NCF_CHECK(nctk_set_collective(ncid,varid))
1321 0 : ABI_MALLOC(num_tasks,(nproc))
1322 0 : do iproc=1,nproc
1323 0 : num_tasks(iproc)=count(mpi_enreg%proc_distrb==iproc-1)
1324 : enddo
1325 0 : num_tasks_max=maxval(num_tasks)
1326 : end if
1327 0 : else if (me==master) then
1328 0 : NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self))
1329 0 : varid=nctk_idname(ncid,"dipole_core_valence")
1330 : !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
1331 : ! NCF_CHECK(nctk_set_collective(ncid,varid))
1332 : !end if
1333 : end if
1334 : end if
1335 :
1336 3 : if (iomode_estf_mpiio) then
1337 : !If MPI-IO, store only elements for one band
1338 12 : ABI_MALLOC(psinablapsi2,(3,nphicor_max,natom,1))
1339 : else
1340 : !If not, store the elements for all bands
1341 0 : ABI_MALLOC(psinablapsi2,(3,nphicor_max,natom,mband))
1342 : end if
1343 15 : pnp_size=size(psinablapsi2)
1344 : !---------------------------------------------------------------------------------
1345 : ! Compute X absorption coefficient and/or X emissivity
1346 :
1347 3 : if (need_absorption) then
1348 15 : ABI_MALLOC(sigx1_k,(nphicor,mom,natom_atnbr))
1349 : end if
1350 3 : if (need_emissivity) then
1351 5 : ABI_MALLOC(emisx_k,(nphicor,mom,natom_atnbr))
1352 : end if
1353 9 : ABI_MALLOC(dhdk2_g,(nphicor))
1354 :
1355 3 : deltae = zero
1356 :
1357 : !LOOP OVER SPINS/K
1358 3 : bdtot_index = 0
1359 3 : itask=0
1360 6 : do isppol=1,nsppol
1361 9 : do ikpt=1,nkpt
1362 3 : etiq=ikpt+(isppol-1)*nkpt
1363 3 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1364 3 : mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))
1365 53 : master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
1366 : ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point
1367 : ! Master node reads and send to relevant processor
1368 3 : if (.not.iomode_estf_mpiio.and.me==master) then
1369 0 : if (iomode==IO_MODE_ETSF) then
1370 0 : nc_start=[1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
1371 0 : nc_count=[3,nphicor_max,natom,mband,1,1]
1372 0 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
1373 : else
1374 0 : psinablapsi2=zero
1375 0 : if (fform2==612) then ! New OPT2 file format
1376 0 : read(opt2_unt) (((psinablapsi2(1,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
1377 0 : read(opt2_unt) (((psinablapsi2(2,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
1378 0 : read(opt2_unt) (((psinablapsi2(3,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom),iband=1,nband_k)
1379 0 : else if (fform2==613) then ! Large OPT2 file format
1380 0 : do iband=1,nband_k
1381 0 : read(opt2_unt) ((psinablapsi2(1,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
1382 0 : read(opt2_unt) ((psinablapsi2(2,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
1383 0 : read(opt2_unt) ((psinablapsi2(3,icor,iatom,iband),icor=1,nphicor_max),iatom=1,natom)
1384 : end do
1385 : else
1386 : !The old writing was not efficient (indexes order is bad)
1387 0 : do iatom=1,natom
1388 0 : read(opt2_unt) ((psinablapsi2(1,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
1389 0 : read(opt2_unt) ((psinablapsi2(2,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
1390 0 : read(opt2_unt) ((psinablapsi2(3,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor_max)
1391 : end do
1392 : end if
1393 : end if
1394 0 : if (.not.mykpt) then
1395 0 : call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,master_band,comm,etiq,ierr)
1396 : end if
1397 : end if
1398 : !! Select k-points for current proc
1399 3 : if (mykpt) then
1400 :
1401 9 : ABI_MALLOC(eig0_k,(nband_k))
1402 6 : ABI_MALLOC(occ_k,(nband_k))
1403 :
1404 41012 : if (need_absorption) sigx1_k=zero
1405 16007 : if (need_emissivity) emisx_k=zero
1406 :
1407 : ! k-dependent data
1408 53 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1409 53 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1410 :
1411 : ! In case of non MPI-IO, receive all (n,m) dipoles from master proc
1412 : ! Then broadcast them to all band processors
1413 3 : if (.not.iomode_estf_mpiio) then
1414 0 : if (me/=master.and.me==master_band) then
1415 0 : call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,me,comm,etiq,ierr)
1416 : end if
1417 0 : call xmpi_bcast(psinablapsi2,master,mpi_enreg%comm_band,mpierr)
1418 : end if
1419 :
1420 : ! LOOP OVER BANDS
1421 :
1422 53 : do iband=1,nband_k
1423 :
1424 : !If MPI-IO, store only ib elements for each iband
1425 : !If not, store all elements
1426 50 : my_iband=merge(1,iband,iomode_estf_mpiio)
1427 240 : dhdk2_g = zero
1428 :
1429 : ! Select bands for current proc
1430 50 : myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me)
1431 53 : if (myband) then
1432 :
1433 : ! In case of MPI-IO, read core-valence dipoles for band n
1434 50 : if (iomode_estf_mpiio) then
1435 50 : itask=itask+1
1436 350 : nc_start=[1,1,1,iband,ikpt,isppol];nc_stride=[1,1,1,1,1,1]
1437 350 : nc_count=[3,nphicor_max,natom,1,1,1]
1438 50 : NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count))
1439 : end if
1440 :
1441 : ! LOOP OVER ATOMS
1442 50 : iatom_atnbr=0
1443 190 : do iatom=1,natom
1444 190 : if(typat(iatom)==itypat_atnbr) then
1445 140 : iatom_atnbr=iatom_atnbr+1
1446 660 : dhdk2_g = zero
1447 660 : do icor=1,nphicor
1448 2220 : do l1=1,3
1449 2080 : dhdk2_g(icor)=dhdk2_g(icor)+psinablapsi2(l1,icor,iatom,my_iband)*third
1450 : end do
1451 : end do
1452 140140 : do iom=1,mom
1453 660140 : do icor=1,nphicor
1454 520000 : diff_occ =occ_cor(icor,itypat_atnbr)/maxocc_cor(icor,itypat_atnbr)-occ_k(iband)/two*nsppol*nspinor
1455 520000 : diff_occ = diff_occ*two/nspinor/nsppol
1456 520000 : diff_eig=eig0_k(iband)-energy_cor(icor,itypat_atnbr)
1457 520000 : oml=oml_edge(icor,iom)
1458 520000 : if(need_absorption) then
1459 520000 : docc_deig=abs(diff_occ/oml)
1460 520000 : if(broad_mode==1) then
1461 520000 : dirac=dom_var1(icor,iom)/((diff_eig-oml)**2+(dom_var1(icor,iom))**2)/pi
1462 : else
1463 0 : dirac=dexp(-((diff_eig-oml)/(sqrt(two)*dom))**2)/(dom*dsqrt(two*pi))
1464 : endif
1465 520000 : if(dirac<1d-20) dirac=zero
1466 520000 : sigx1_k(icor,iom,iatom_atnbr)=sigx1_k(icor,iom,iatom_atnbr)+dhdk2_g(icor)*docc_deig*dirac*pi/ucvol
1467 : endif
1468 660000 : if (need_emissivity) then
1469 120000 : docc_deig=abs(occ_k(iband)/oml)
1470 120000 : if(broad_mode==1) then
1471 120000 : dirac=dom/((diff_eig-oml)**2+dom**2)/pi
1472 : else
1473 0 : dirac=dexp(-((diff_eig-oml)/(dom*sqrt(two)))**2)/(dom*dsqrt(two*pi))
1474 : endif
1475 120000 : if(dirac<1d-20) dirac=zero
1476 120000 : emisx_k(icor,iom,iatom_atnbr)=emisx_k(icor,iom,iatom_atnbr)+dhdk2_g(icor)*docc_deig*dirac*pi/ucvol
1477 : endif
1478 : end do
1479 : end do
1480 : endif
1481 : end do ! iatom
1482 : end if ! my band?
1483 : end do ! iband
1484 :
1485 : ! Accumulate k-point contribution
1486 3 : if (need_absorption) then
1487 : sigx1(1:nphicor,1:mom,1:natom_atnbr,isppol)=sigx1(1:nphicor,1:mom,1:natom_atnbr,isppol) &
1488 41012 : & +wtk(ikpt)*sigx1_k(1:nphicor,1:mom,1:natom_atnbr)
1489 : end if
1490 3 : if (need_emissivity) then
1491 : emisx(1:nphicor,1:mom,1:natom_atnbr,isppol)=emisx(1:nphicor,1:mom,1:natom_atnbr,isppol) &
1492 16005 : & +wtk(ikpt)*emisx_k(1:nphicor,1:mom,1:natom_atnbr)
1493 : end if
1494 :
1495 : ! Validity limit
1496 3 : deltae=deltae+eig0_k(nband_k)
1497 :
1498 3 : ABI_FREE(eig0_k)
1499 3 : ABI_FREE(occ_k)
1500 : ! End loop over kpt/spin
1501 : end if ! My kpt?
1502 6 : bdtot_index=bdtot_index+nband_k
1503 : end do ! ikpt
1504 : end do ! isppol
1505 :
1506 3 : if (iomode == IO_MODE_ETSF) then
1507 3 : if(iomode_estf_mpiio.and.nproc>1) then
1508 0 : do idum=num_tasks(me+1)+1,num_tasks_max
1509 0 : nc_start=[1,1,1,1,1,1];nc_stride=[1,1,1,1,1,1]
1510 0 : nc_count=[0,0,0,0,0,0]
1511 0 : NCF_CHECK(nf90_get_var(ncid,varid,dummy,start=nc_start,stride=nc_stride,count=nc_count))
1512 : enddo
1513 0 : ABI_FREE(num_tasks)
1514 : endif
1515 : endif
1516 :
1517 : !Accumulate kpt/band contributions over processors
1518 3 : if (need_absorption) then
1519 3 : call xmpi_sum(sigx1,comm,mpierr)
1520 : end if
1521 3 : if (need_emissivity) then
1522 1 : call xmpi_sum(emisx,comm,mpierr)
1523 : end if
1524 3 : call xmpi_sum(deltae,comm,mpierr)
1525 3 : deltae=deltae/mpi_enreg%nproc_band
1526 :
1527 :
1528 : !Release some memory
1529 3 : ABI_FREE(dhdk2_g)
1530 3 : ABI_FREE(psinablapsi2)
1531 3 : if (need_absorption) then
1532 3 : ABI_FREE(sigx1_k)
1533 : end if
1534 3 : if (need_emissivity) then
1535 1 : ABI_FREE(emisx_k)
1536 : end if
1537 :
1538 : !Close core-valence dipoles file
1539 3 : if (iomode == IO_MODE_ETSF) then
1540 3 : if (iomode_estf_mpiio.or.me==master) then
1541 3 : NCF_CHECK(nf90_close(ncid))
1542 : end if
1543 0 : else if (me==master) then
1544 0 : ierr=close_unit(opt2_unt,msg)
1545 0 : ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam2))
1546 : end if
1547 :
1548 : !---------------------------------------------------------------------------------
1549 : ! Post-processing
1550 :
1551 3 : if (need_absorption) then
1552 41015 : sigx1=sigx1*dble(natom_atnbr)
1553 15 : ABI_MALLOC(sigx1_av,(nphicor,mom,nsppol))
1554 14006 : sigx1_av=zero
1555 : endif
1556 3 : if (need_emissivity) then
1557 : ! Filter low values
1558 2 : do isppol=1,nsppol
1559 6 : do iatom_atnbr=1,natom_atnbr
1560 4005 : do iom=1,mom
1561 16004 : do icor=1,nphicor
1562 16000 : if (emisx(icor,iom,iatom_atnbr,isppol)<=tol16) emisx(icor,iom,iatom_atnbr,isppol)=zero
1563 : end do
1564 : end do
1565 : end do
1566 : end do
1567 16006 : emisx=emisx*dble(natom_atnbr)
1568 5 : ABI_MALLOC(emisx_av,(nphicor,mom,nsppol))
1569 4002 : emisx_av=zero
1570 : endif
1571 :
1572 6 : do isppol=1,nsppol
1573 15 : do iatom_atnbr=1,natom_atnbr
1574 9012 : do iom=1,mom
1575 41009 : do icor=1,nphicor
1576 32000 : if(need_absorption) sigx1_av(icor,iom,isppol)=sigx1_av(icor,iom,isppol)+sigx1(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
1577 41000 : if(need_emissivity) emisx_av(icor,iom,isppol)=emisx_av(icor,iom,isppol)+emisx(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
1578 : end do
1579 : end do
1580 : end do
1581 : end do
1582 :
1583 : ! Spin treatment
1584 3 : if(nsppol==2) then
1585 0 : if(need_absorption) then
1586 0 : ABI_MALLOC(sum_spin_sigx1,(nphicor,mom,natom_atnbr))
1587 0 : ABI_MALLOC(sum_spin_sigx1_av,(nphicor,mom))
1588 0 : sum_spin_sigx1=zero ; sum_spin_sigx1_av=zero
1589 : endif
1590 0 : if(need_emissivity) then
1591 0 : ABI_MALLOC(sum_spin_emisx,(nphicor,mom,natom_atnbr))
1592 0 : ABI_MALLOC(sum_spin_emisx_av,(nphicor,mom))
1593 0 : sum_spin_emisx=zero ; sum_spin_emisx_av=zero
1594 : endif
1595 0 : do isppol=1,nsppol
1596 0 : do iatom_atnbr=1,natom_atnbr
1597 0 : do iom=1,mom
1598 0 : do icor=1,nphicor
1599 0 : if(need_absorption) sum_spin_sigx1(icor,iom,iatom_atnbr)=sum_spin_sigx1(icor,iom,iatom_atnbr) &
1600 0 : & +sigx1(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
1601 0 : if(need_emissivity) sum_spin_emisx(icor,iom,iatom_atnbr)=sum_spin_emisx(icor,iom,iatom_atnbr) &
1602 0 : & +emisx(icor,iom,iatom_atnbr,isppol)/dble(natom_atnbr)
1603 : end do
1604 : end do
1605 : end do
1606 : end do
1607 0 : do isppol=1,nsppol
1608 0 : do iom=1,mom
1609 0 : do icor=1,nphicor
1610 0 : if(need_absorption) sum_spin_sigx1_av(icor,iom)=sum_spin_sigx1_av(icor,iom)+sigx1_av(icor,iom,isppol)
1611 0 : if(need_emissivity) sum_spin_emisx_av(icor,iom)=sum_spin_emisx_av(icor,iom)+emisx_av(icor,iom,isppol)
1612 : end do
1613 : end do
1614 : end do
1615 : endif
1616 :
1617 : ! Units
1618 3 : if(au_units==0) then
1619 14003 : oml_edge=oml_edge*Ha_eV
1620 3 : if(need_absorption) then
1621 14006 : sigx1_av=sigx1_av*Ohmcm
1622 41015 : sigx1=sigx1*Ohmcm
1623 3 : if(nsppol==2) then
1624 0 : sum_spin_sigx1_av=sum_spin_sigx1_av*Ohmcm
1625 0 : sum_spin_sigx1=sum_spin_sigx1*Ohmcm
1626 : endif
1627 : endif
1628 3 : if(need_emissivity) then
1629 4002 : emisx_av=emisx_av*Ohmcm
1630 16006 : emisx=emisx*Ohmcm
1631 1 : if(nsppol==2) then
1632 0 : sum_spin_emisx_av=sum_spin_emisx_av*Ohmcm
1633 0 : sum_spin_emisx=sum_spin_emisx*Ohmcm
1634 : endif
1635 : endif
1636 : endif
1637 :
1638 : !---------------------------------------------------------------------------------
1639 : ! Output results
1640 3 : write(str_atm,*) input_atm
1641 3 : str_atm=adjustl(str_atm)
1642 : str_atm=trim(str_atm)
1643 :
1644 : !Only master node outputs results in files (only master node)
1645 3 : if (me==master) then
1646 :
1647 : ! Standard output
1648 3 : if (need_absorption) then
1649 3 : write(std_out,*) 'Absorption: valence state orbital energies: omin,omax',omin_sig,omax_sig
1650 : end if
1651 3 : if (need_emissivity) then
1652 1 : write(std_out,*) 'Emissivity: valence state orbital energies: omin,omax',omin_emis,omax_emis
1653 : end if
1654 3 : if (need_absorption) then
1655 : write(std_out,'(a,f10.5,a,f10.5,a)')&
1656 3 : & ' Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
1657 : end if
1658 :
1659 : ! _sigX file
1660 : if (need_absorption) then
1661 3 : if (open_file(trim(filnam_out)//'_sigX_at'//str_atm,msg,newunit=sigx1_unt,form='formatted',action="write")/=0) then
1662 0 : ABI_ERROR(msg)
1663 : end if
1664 3 : if (abs(dom_max)>tol10) then
1665 1 : write(sigx1_unt,'(a)')'#***************************************************** ARCTAN SMEARING ********'
1666 1 : write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1667 1 : write(sigx1_unt,'(a,2f10.5,a)')'# dom_max,center =',dom_max,dom_ctr,' Ha'
1668 : else
1669 2 : write(sigx1_unt,'(a)')'#***************************************************** FIXED SMEARING ********'
1670 2 : write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha'
1671 : endif
1672 3 : write(sigx1_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1)
1673 3 : write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2)
1674 3 : write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3)
1675 3 : write(sigx1_unt,'(a,i8)' ) '# natom =',natom
1676 3 : write(sigx1_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol
1677 3 : write(sigx1_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha'
1678 3 : write(sigx1_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV'
1679 3 : write(sigx1_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin'
1680 3 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
1681 3 : write(sigx1_unt,'(a,i4)') '# Number of core orbitals nc=',nphicor
1682 14 : do icor=1,nphicor
1683 14 : if (kappacor(icor,itypat_atnbr)==0) then
1684 3 : write(sigx1_unt,'(a,2i4,4f15.5)') '# n, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ', &
1685 3 : ncor(icor,itypat_atnbr),lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
1686 6 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1687 : else
1688 8 : if (kappacor(icor,itypat_atnbr)>0) then
1689 2 : j2=2*lcor(icor,itypat_atnbr)-1
1690 2 : write(sigx1_unt,'(a,i4,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1691 2 : & ncor(icor,itypat_atnbr),j2,' /2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
1692 4 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1693 : else
1694 6 : if (kappacor(icor,itypat_atnbr)<-1) then
1695 2 : j2=2*lcor(icor,itypat_atnbr)+1
1696 2 : write(sigx1_unt,'(a,i4,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1697 2 : & ncor(icor,itypat_atnbr),j2,'/2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
1698 4 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1699 : else
1700 4 : write(sigx1_unt,'(a,i4,a,i4,4f15.5)') '# n, j, l, Energy(Ha), Edge(Ha), Energy(eV), Edge(eV): ',&
1701 4 : & ncor(icor,itypat_atnbr),'1/2',lcor(icor,itypat_atnbr),energy_cor(icor,itypat_atnbr),edge(icor),&
1702 8 : & energy_cor(icor,itypat_atnbr)*Ha_eV,edge(icor)*Ha_eV
1703 : end if
1704 : end if
1705 : end if
1706 : end do
1707 3 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
1708 : write(sigx1_unt,'(a,f10.5,a,f10.5,a)')&
1709 3 : & '# Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
1710 3 : write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------'
1711 3 : if(au_units>0) then
1712 0 : write(sigx1_unt,'(a)')'# om(au),sig1_av(au),sig1(au) for each orbital'
1713 : else
1714 3 : write(sigx1_unt,'(a)')'# om(eV),sig1_av(Ohm.cm)-1,sig1(Ohm.cm)-1 for each orbital'
1715 : endif
1716 3003 : do iom=1,mom
1717 3003 : if(nsppol==1) then
1718 : write(sigx1_unt,'(100(1x,e15.8))') &
1719 14000 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1),sigx1(icor,iom,atnbr,1),icor=1,nphicor)
1720 : else
1721 : write(sigx1_unt,'(100(1x,e15.8))') &
1722 0 : & (oml_edge(icor,iom),sum_spin_sigx1_av(icor,iom),sum_spin_sigx1(icor,iom,atnbr),icor=1,nphicor)
1723 : endif
1724 : end do
1725 3 : close(sigx1_unt)
1726 : ! absX file
1727 3 : if(au_units==0) then
1728 3 : if(open_file(trim(filnam_out)//'_absX_at'//str_atm,msg,newunit=absx_unt,form='formatted',action="write")/=0) then
1729 0 : ABI_ERROR(msg)
1730 : end if
1731 3 : write(absx_unt,'(a)') '# om(eV), abso X average (cm-1), abso X(cm-1)'
1732 3003 : do iom=1,mom
1733 3003 : if(nsppol==1) then
1734 : write(absx_unt,'(100(1x,e15.8))') &
1735 11000 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1)/(Speed_Light_SI*eps0),&
1736 14000 : & sigx1(icor,iom,atnbr,1)/(Speed_Light_SI*eps0),icor=1,nphicor)
1737 : else
1738 : write(absx_unt,'(100(1x,e15.8))') &
1739 0 : & (oml_edge(icor,iom),sum_spin_sigx1_av(icor,iom)/(Speed_Light_SI*eps0),&
1740 0 : & sum_spin_sigx1(icor,iom,atnbr)/(Speed_Light_SI*eps0),icor=1,nphicor)
1741 : endif
1742 : end do
1743 3 : close(absx_unt)
1744 : endif
1745 : end if
1746 :
1747 : ! _s_sigX
1748 3 : if (need_absorption.and.nsppol==2) then
1749 0 : if (open_file(trim(filnam_out)//'_sigX_up_at'//str_atm,msg,newunit=sigx1_up_unt,form='formatted',action="write")/=0) then
1750 0 : ABI_ERROR(msg)
1751 : end if
1752 0 : if (open_file(trim(filnam_out)//'_sigX_dn_at'//str_atm,msg,newunit=sigx1_dn_unt,form='formatted',action="write")/=0) then
1753 0 : ABI_ERROR(msg)
1754 : end if
1755 0 : do iom=1,mom
1756 : write(sigx1_up_unt,'(100(1x,e15.8))') &
1757 0 : & (oml_edge(icor,iom),sigx1_av(icor,iom,1),sigx1(icor,iom,atnbr,1),icor=1,nphicor)
1758 : write(sigx1_dn_unt,'(100(1x,e15.8))') &
1759 0 : & (oml_edge(icor,iom),sigx1_av(icor,iom,2),sigx1(icor,iom,atnbr,2),icor=1,nphicor)
1760 : end do
1761 0 : close(sigx1_up_unt)
1762 0 : close(sigx1_dn_unt)
1763 : end if
1764 :
1765 : ! _emisX file
1766 3 : if (need_emissivity) then
1767 1 : if (open_file(trim(filnam_out)//'_emisX_at'//str_atm,msg,newunit=ems_unt,form='formatted',action="write")/=0) then
1768 0 : ABI_ERROR(msg)
1769 : end if
1770 1 : if(au_units>0) then
1771 0 : write(ems_unt,'(a)')'# om(au),sig1_av(au),sig1(au) for each orbital'
1772 : else
1773 1 : write(ems_unt,'(a)')'# om(eV),sig1_av(Ohm.cm)-1,sig1(Ohm.cm)-1 for each orbital'
1774 : endif
1775 1001 : do iom=1,mom
1776 1001 : if(nsppol==1) then
1777 : write(ems_unt,'(3(3(1x,e15.8),2x))') &
1778 4000 : & (oml_edge(icor,iom),emisx_av(icor,iom,1),emisx(icor,iom,atnbr,1),icor=1,nphicor)
1779 : else
1780 : write(ems_unt,'(3(3(1x,e15.8),2x))') &
1781 0 : & (oml_edge(icor,iom),sum_spin_emisx_av(icor,iom)/dble(natom_atnbr),sum_spin_emisx(icor,iom,atnbr),icor=1,nphicor)
1782 : endif
1783 : end do
1784 1 : close(ems_unt)
1785 : end if
1786 :
1787 : ! _s_emisX
1788 1 : if (need_emissivity.and.nsppol==2) then
1789 0 : if (open_file(trim(filnam_out)//'_emisX_up_at'//str_atm,msg,newunit=ems_up_unt,form='formatted',action="write")/=0) then
1790 0 : ABI_ERROR(msg)
1791 : end if
1792 0 : if (open_file(trim(filnam_out)//'_emisX_dn_at'//str_atm,msg,newunit=ems_dn_unt,form='formatted',action="write")/=0) then
1793 0 : ABI_ERROR(msg)
1794 : end if
1795 0 : do iom=1,mom
1796 : write(ems_up_unt,'(3(3(1x,e15.8),2x))') &
1797 0 : & (oml_edge(icor,iom),emisx_av(icor,iom,1)/dble(natom_atnbr),emisx(icor,iom,atnbr,1),icor=1,nphicor)
1798 : write(ems_dn_unt,'(3(3(1x,e15.8),2x))') &
1799 0 : (oml_edge(icor,iom),emisx_av(icor,iom,2)/dble(natom_atnbr),emisx(icor,iom,atnbr,2),icor=1,nphicor)
1800 : end do
1801 0 : close(ems_up_unt)
1802 0 : close(ems_dn_unt)
1803 : end if
1804 :
1805 : endif ! master node
1806 :
1807 : !---------------------------------------------------------------------------------
1808 : ! End
1809 :
1810 : !Release memory space
1811 3 : if (need_absorption) then
1812 3 : ABI_FREE(sigx1)
1813 3 : ABI_FREE(sigx1_av)
1814 3 : if (nsppol==2) then
1815 0 : ABI_FREE(sum_spin_sigx1)
1816 0 : ABI_FREE(sum_spin_sigx1_av)
1817 : end if
1818 3 : ABI_FREE(dom_var1)
1819 3 : ABI_FREE(oml_edge)
1820 : end if
1821 3 : if (need_emissivity) then
1822 1 : ABI_FREE(emisx)
1823 1 : ABI_FREE(emisx_av)
1824 1 : if (nsppol==2) then
1825 0 : ABI_FREE(sum_spin_emisx)
1826 0 : ABI_FREE(sum_spin_emisx_av)
1827 : end if
1828 1 : ABI_FREE(oml_emis)
1829 : end if
1830 3 : ABI_FREE(typat)
1831 3 : ABI_FREE(ncor)
1832 3 : ABI_FREE(lcor)
1833 3 : ABI_FREE(kappacor)
1834 3 : ABI_FREE(energy_cor)
1835 3 : ABI_FREE(edge)
1836 3 : ABI_FREE(eigen0)
1837 3 : ABI_FREE(nband)
1838 3 : ABI_FREE(occ)
1839 3 : ABI_FREE(occ_cor)
1840 3 : ABI_FREE(maxocc_cor)
1841 3 : ABI_FREE(wtk)
1842 3 : call hdr%free()
1843 3 : call destroy_mpi_enreg(mpi_enreg)
1844 :
1845 6 : end subroutine conducti_paw_core
1846 : !!***
1847 :
1848 : !----------------------------------------------------------------------
1849 :
1850 : !!****f* m_conducti/conducti_nc
1851 : !! NAME
1852 : !! conducti_nc
1853 : !!
1854 : !! FUNCTION
1855 : !! This program computes the elements of the optical frequency dependent
1856 : !! conductivity tensor and the conductivity along the three principal axes
1857 : !! from the Kubo-Greenwood formula.
1858 : !!
1859 : !! INPUTS
1860 : !! (main routine)
1861 : !!
1862 : !! OUTPUT
1863 : !! (main routine)
1864 : !!
1865 : !! NOTES
1866 : !! bantot
1867 : !! doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy.
1868 : !! dom=frequency range
1869 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
1870 : !! eigen11(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
1871 : !! in reciprocal direction 100
1872 : !! eigen12(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
1873 : !! in reciprocal direction 010
1874 : !! eigen13(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree)
1875 : !! in reciprocal direction 001
1876 : !! ecut=kinetic energy planewave cutoff (hartree).
1877 : !! entropy= entropy associated with the smearing (adimensional)
1878 : !! fermie= fermi energy (Hartree)
1879 : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$).
1880 : !! gmet_inv(3,3)=inverse of reciprocal space metric ($\textrm{bohr}^{2}$).
1881 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
1882 : !! kin11= Onsager kinetic coeficient=optical conductivity
1883 : !! kin12= Onsager kinetic coeficient
1884 : !! kin21= Onsager kinetic coeficient
1885 : !! kin22= Onsager kinetic coeficient
1886 : !! Kth=thermal conductivity
1887 : !! mom=number of frequency for conductivity computation
1888 : !! mband=maximum number of bands.
1889 : !! natom = number of atoms in the unit cell.
1890 : !! nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
1891 : !! nelect=number of electrons per unit cell
1892 : !! nkpt=number of k points in the IBZ for this perturbation
1893 : !! ngfft(3)=integer fft box dimensions.
1894 : !! nspinor=number of spinorial components of the wavefunctions.
1895 : !! nsppol=1 for unpolarized, 2 for spin-polarized.
1896 : !! ntypat = number of atom types.
1897 : !! occ(mband*nkpt*nsppol)=occupation number for each band and k.
1898 : !! occopt==option for occupancies
1899 : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
1900 : !! rprimd(3,3)=real space primitive translations.
1901 : !! of primitive translations.
1902 : !! Sth=thermopower
1903 : !! tsmear=smearing width (or temperature) in Hartree
1904 : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$).
1905 : !! wind=frequency windows for computations of sigma
1906 : !! wtk(nkpt)=weight assigned to each k point.
1907 : !! znucl(natom)=atomic number of atoms
1908 : !! np_sum=noziere-pines sumrule
1909 : !! cond_kg(mom)=kubo-greenwood conductivity
1910 : !!
1911 : !! SOURCE
1912 :
1913 :
1914 1 : subroutine conducti_nc(filnam,filnam_out)
1915 :
1916 : !Arguments -----------------------------------
1917 : !scalars
1918 : character(len=fnlen),intent(in) :: filnam,filnam_out
1919 :
1920 : !Local variables-------------------------------
1921 : !scalars
1922 : integer,parameter :: formeig0=0,formeig1=1
1923 : integer :: bantot,bd2tot_index,bdtot0_index,bdtot_index
1924 : integer :: headform,iband,ii,jj,ikpt,iunt
1925 : integer :: index_1,iom,isppol,jband,l1,l2,mband,mom,natom,nband1
1926 : integer :: nrot,iomode
1927 : integer :: nband_k,nkpt,nlign,nrest,nspinor,nsppol,ntypat
1928 : integer :: occopt,comm
1929 : integer :: tens_unt,lij_unt,sig_unt,kth_unt,ocond_unt
1930 : real(dp) :: deltae,dosdeltae,diff_occ,dom,ecut,entropy,fermie,maxocc
1931 : real(dp) :: nelect,np_sum,np_sum_k1,np_sum_k2,omin,oml,socc,socc_k,sig
1932 : real(dp) :: tphysel,tsmear,ucvol,wind,Tatm
1933 : character(len=fnlen) :: filnam0,filnam1,filnam2,filnam3
1934 : character(len=500) :: msg
1935 1 : type(hdr_type) :: hdr
1936 14 : type(wfk_t) :: gswfk,ddk1,ddk2,ddk3
1937 : !arrays
1938 1 : integer,allocatable :: nband(:)
1939 : real(dp) :: gmet(3,3),gmet_inv(3,3),gprimd(3,3),gprimd_inv(3,3),rmet(3,3),rprimd(3,3)
1940 2 : real(dp),allocatable :: cond_kg(:,:,:),cond_kg_cart(:,:,:),cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:)
1941 3 : real(dp),allocatable :: doccde(:),doccde_k(:),cond_kg_xx(:),cond_kg_yy(:),cond_kg_zz(:),trace(:)
1942 1 : real(dp),allocatable :: eig0_k(:),eig0tmp(:),eig1_k(:,:),eigen0(:),eigen11(:)
1943 1 : real(dp),allocatable :: eigen12(:),eigtmp(:)
1944 2 : real(dp),allocatable :: eigen13(:),occ(:),occ_k(:),wtk(:),cond_tot(:),oml1(:)
1945 1 : real(dp),allocatable :: kin11(:),kin12(:),kin21(:),kin22(:)
1946 1 : real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
1947 : real(dp) :: cond_kg_w(3,3),z(3,3)
1948 : real(dp) :: eig_cond(3)
1949 :
1950 : ! *********************************************************************************
1951 :
1952 : !Read data file
1953 1 : if (open_file(filnam,msg,newunit=iunt,form='formatted',status="old")/=0) then
1954 0 : ABI_ERROR(msg)
1955 : end if
1956 :
1957 1 : rewind(iunt)
1958 1 : read(iunt,*)
1959 1 : read(iunt,'(a)')filnam1 ! first ddk file
1960 1 : read(iunt,'(a)')filnam2 ! second ddk file
1961 1 : read(iunt,'(a)')filnam3 ! third ddk file
1962 1 : read(iunt,'(a)')filnam0 ! ground-state data
1963 :
1964 : !Open the GS Wavefunction file and the 3 DDK files.
1965 : ! TODO: one should perform basic consistency tests for the GS WFK and the DDK files, e.g.
1966 : ! k-points and their order, spins, number of bands could differ in the four files.
1967 : ! Note indeed that we are assuming the same numer of bands in all the files.
1968 1 : comm = xmpi_comm_self
1969 1 : call nctk_fort_or_ncfile(filnam0, iomode, msg)
1970 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
1971 1 : call gswfk%open_read(filnam0, formeig0, iomode, get_unit(), comm)
1972 :
1973 1 : call nctk_fort_or_ncfile(filnam1, iomode, msg)
1974 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
1975 1 : call ddk1%open_read(filnam1, formeig1, iomode, get_unit(), comm, hdr_out=hdr)
1976 :
1977 1 : call nctk_fort_or_ncfile(filnam2, iomode, msg)
1978 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
1979 1 : call ddk2%open_read(filnam2, formeig1, iomode, get_unit(), comm)
1980 :
1981 1 : call nctk_fort_or_ncfile(filnam3, iomode, msg)
1982 1 : if (len_trim(msg) /= 0) ABI_ERROR(msg)
1983 1 : call ddk3%open_read(filnam3, formeig1, iomode, get_unit(), comm)
1984 :
1985 1 : if (ddk1%compare(ddk2) /= 0) then
1986 0 : ABI_ERROR("ddk1 and ddk2 are not consistent. see above messages")
1987 : end if
1988 1 : if (ddk1%compare(ddk3) /= 0) then
1989 0 : ABI_ERROR("ddk1 and ddk3 are not consistent. see above messages")
1990 : end if
1991 :
1992 : !Extract params from the header of the first ddk file (might have been the GS file ?)
1993 :
1994 : !Extract info from the header
1995 1 : headform=hdr%headform
1996 1 : bantot=hdr%bantot
1997 1 : ecut=hdr%ecut_eff
1998 1 : natom=hdr%natom
1999 1 : nkpt=hdr%nkpt
2000 1 : nspinor=hdr%nspinor
2001 1 : nsppol=hdr%nsppol
2002 1 : ntypat=hdr%ntypat
2003 1 : occopt=hdr%occopt
2004 13 : rprimd(:,:)=hdr%rprimd(:,:)
2005 3 : ABI_MALLOC(nband,(nkpt*nsppol))
2006 3 : ABI_MALLOC(occ,(bantot))
2007 1 : fermie=hdr%fermie
2008 21 : occ(1:bantot)=hdr%occ(1:bantot)
2009 2 : nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
2010 :
2011 : !Get mband, as the maximum value of nband(nkpt)
2012 2 : mband=maxval(nband(:))
2013 :
2014 1 : write(std_out,*)
2015 1 : write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1)
2016 1 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2)
2017 1 : write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3)
2018 1 : write(std_out,'(a,i8)') ' natom =',natom
2019 1 : write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband
2020 1 : write(std_out,'(a, f10.5,a)' ) ' ecut =',ecut,' Ha'
2021 1 : write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV'
2022 :
2023 : !Prepare the reading of ddk Wff files
2024 3 : ABI_MALLOC(eigtmp,(2*mband*mband))
2025 3 : ABI_MALLOC(eig0tmp,(mband))
2026 :
2027 : !Read the eigenvalues of ground-state and ddk files
2028 3 : ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
2029 3 : ABI_MALLOC(eigen11,(2*mband*mband*nkpt*nsppol))
2030 2 : ABI_MALLOC(eigen12,(2*mband*mband*nkpt*nsppol))
2031 2 : ABI_MALLOC(eigen13,(2*mband*mband*nkpt*nsppol))
2032 1 : bdtot0_index=0 ; bdtot_index=0
2033 2 : do isppol=1,nsppol
2034 3 : do ikpt=1,nkpt
2035 1 : nband1=nband(ikpt+(isppol-1)*nkpt)
2036 1 : call gswfk%read_eigk(ikpt,isppol,xmpio_single,eig0tmp)
2037 21 : eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1)
2038 :
2039 1 : call ddk1%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2040 801 : eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2041 :
2042 1 : call ddk2%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2043 801 : eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2044 :
2045 1 : call ddk3%read_eigk(ikpt,isppol,xmpio_single,eigtmp)
2046 801 : eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
2047 :
2048 1 : bdtot0_index=bdtot0_index+nband1
2049 2 : bdtot_index=bdtot_index+2*nband1**2
2050 : end do
2051 : end do
2052 :
2053 : !Close files
2054 1 : call gswfk%close()
2055 1 : call ddk1%close()
2056 1 : call ddk2%close()
2057 1 : call ddk3%close()
2058 :
2059 1 : ABI_FREE(eigtmp)
2060 1 : ABI_FREE(eig0tmp)
2061 :
2062 : !---------------------------------------------------------------------------------
2063 : !Gmet inversion
2064 1 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2065 1 : call matr3inv(gmet,gmet_inv)
2066 1 : call matr3inv(gprimd,gprimd_inv)
2067 :
2068 : !---------------------------------------------------------------------------------
2069 : !Derivative of occupation wrt the energy.
2070 :
2071 2 : ABI_MALLOC(doccde,(mband*nkpt*nsppol))
2072 3 : ABI_MALLOC(wtk,(nkpt))
2073 :
2074 1 : read(iunt,*)tsmear
2075 1 : Tatm=tsmear*Ha_K
2076 1 : write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin'
2077 : !
2078 1 : nlign=nkpt/6
2079 1 : nrest=nkpt-6*nlign
2080 1 : index_1=0
2081 1 : do ii=1,nlign
2082 0 : read(iunt,*)wtk(1+index_1:6+index_1)
2083 1 : index_1=index_1+6
2084 : end do
2085 1 : if (nrest/=0) then
2086 1 : read(iunt,*)wtk(6*nlign+1:nkpt)
2087 : end if
2088 : !
2089 1 : if (occopt==1) then
2090 0 : write(std_out,'(a,i4)') ' occopt =',occopt
2091 0 : doccde=zero
2092 : else
2093 1 : tphysel=zero
2094 1 : maxocc=two/(nsppol*nspinor)
2095 1 : dosdeltae=zero
2096 : ! CP: using 1 and nband(0) as dummy value, because function
2097 : ! not implemented for occopt==9; adding fermih=fermie in the list of arguments as well
2098 : call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermie,maxocc,mband,nband,&
2099 1 : & nelect,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,11,wtk,1,nband(1))
2100 : end if
2101 :
2102 : !---------------------------------------------------------------------------------
2103 : !Size of the frequency range
2104 :
2105 1 : read(iunt,*)dom,wind
2106 1 : close(iunt)
2107 1 : mom=int(wind/dom)
2108 3 : ABI_MALLOC(oml1,(mom))
2109 273 : do iom=1,mom
2110 273 : oml1(iom)=tol10*1000._dp+dble(iom)*dom
2111 : end do
2112 :
2113 4 : ABI_MALLOC(cond_nd,(mom,3,3))
2114 2 : ABI_MALLOC(cond_kg,(mom,3,3))
2115 2 : ABI_MALLOC(cond_kg_cart,(mom,3,3))
2116 2 : ABI_MALLOC(cond_kg_xx,(mom))
2117 2 : ABI_MALLOC(cond_kg_yy,(mom))
2118 2 : ABI_MALLOC(trace,(mom))
2119 2 : ABI_MALLOC(cond_kg_zz,(mom))
2120 2 : ABI_MALLOC(cond_tot,(mom))
2121 2 : ABI_MALLOC(kin11,(mom))
2122 2 : ABI_MALLOC(kin12,(mom))
2123 2 : ABI_MALLOC(kin21,(mom))
2124 2 : ABI_MALLOC(kin22,(mom))
2125 2 : ABI_MALLOC(kin11_k,(mom))
2126 2 : ABI_MALLOC(kin12_k,(mom))
2127 2 : ABI_MALLOC(kin21_k,(mom))
2128 2 : ABI_MALLOC(kin22_k,(mom))
2129 2 : ABI_MALLOC(Kth,(mom))
2130 2 : ABI_MALLOC(Stp,(mom))
2131 1 : write(std_out,'(a,i8,2f10.5,a)')' mom,wind,dom =',mom,wind,dom,' Ha'
2132 :
2133 : !---------------------------------------------------------------------------------
2134 :
2135 273 : kin11 = zero
2136 273 : kin12 = zero
2137 273 : kin21 = zero
2138 273 : kin22 = zero
2139 1 : np_sum = zero
2140 1 : socc = zero
2141 2461 : cond_kg = zero
2142 :
2143 : !LOOP OVER SPINS
2144 2 : do isppol=1,nsppol
2145 : !
2146 1 : bdtot_index = 0
2147 1 : bd2tot_index = 0
2148 :
2149 1 : deltae = zero
2150 : !
2151 : ! BIG FAT k POINT LOOP
2152 : !
2153 2 : do ikpt=1,nkpt
2154 :
2155 1 : nband_k=nband(ikpt+(isppol-1)*nkpt)
2156 :
2157 3 : ABI_MALLOC(eig0_k,(nband_k))
2158 4 : ABI_MALLOC(eig1_k,(2*nband_k**2,3))
2159 2 : ABI_MALLOC(occ_k,(nband_k))
2160 2 : ABI_MALLOC(doccde_k,(nband_k))
2161 4 : ABI_MALLOC(dhdk2_r,(3,3,nband_k,nband_k))
2162 4 : ABI_MALLOC(dhdk2_g,(nband_k,nband_k))
2163 :
2164 2461 : cond_nd = zero
2165 273 : kin11_k = zero
2166 273 : kin12_k = zero
2167 273 : kin21_k = zero
2168 273 : kin22_k = zero
2169 21 : np_sum_k1 = zero
2170 21 : np_sum_k2 = zero
2171 21 : socc_k = zero
2172 5221 : dhdk2_r = zero
2173 421 : dhdk2_g = zero
2174 :
2175 : ! eigenvalue for k-point
2176 21 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
2177 : ! first derivative eigenvalues for k-point
2178 801 : eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2179 801 : eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2180 801 : eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_k**2+bd2tot_index)
2181 : ! occupation numbers for k-point
2182 21 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
2183 : ! derivative of occupation number for k-point
2184 21 : doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
2185 :
2186 : ! LOOP OVER BAND
2187 21 : do iband=1,nband_k
2188 420 : do jband=1,nband_k
2189 : !
2190 : ! TODO : replace with BLAS calls
2191 1600 : do l1=1,3
2192 5200 : do l2=1,3
2193 15600 : do ii=1,3
2194 46800 : do jj=1,3
2195 : dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)+(rprimd(l1,ii)&
2196 : & *eig1_k(2*iband-1+(jband-1)*2*nband_k,ii)*&
2197 : & rprimd(l2,jj)*eig1_k(2*iband-1+(jband-1)*2*nband_k,jj)&
2198 : & +rprimd(l1,ii)*eig1_k(2*iband +(jband-1)*2*nband_k,ii)*&
2199 43200 : & rprimd(l2,jj)*eig1_k(2*iband+(jband-1)*2*nband_k,jj))
2200 : end do
2201 : end do
2202 : end do
2203 : end do
2204 :
2205 1600 : do l1=1,3
2206 5200 : do l2=1,3
2207 4800 : dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)/two_pi/two_pi
2208 : end do
2209 : end do
2210 :
2211 : ! TODO: replace with BLAS calls
2212 1600 : do l1=1,3
2213 5200 : do l2=1,3
2214 : dhdk2_g(iband,jband)=dhdk2_g(iband,jband)+gmet_inv(l1,l2)*( &
2215 : & eig1_k(2*iband-1+(jband-1)*2*nband_k,l1)*&
2216 : & eig1_k(2*iband-1+(jband-1)*2*nband_k,l2) &
2217 : & +eig1_k(2*iband +(jband-1)*2*nband_k,l1)*&
2218 4800 : & eig1_k(2*iband +(jband-1)*2*nband_k,l2))
2219 : end do
2220 : end do
2221 400 : dhdk2_g(iband,jband)=dhdk2_g(iband,jband)/two_pi/two_pi
2222 :
2223 400 : diff_occ = occ_k(iband)-occ_k(jband)
2224 : ! if (dabs(diff_occ)>=tol8) then
2225 :
2226 : ! Conductivity for each omega
2227 400 : omin = zero
2228 109200 : do iom=1,mom
2229 108800 : oml=oml1(iom)
2230 108800 : if (jband>iband) then
2231 : sig= dhdk2_g(iband,jband)&
2232 : & *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)&
2233 51680 : & -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2))
2234 51680 : kin11_k(iom)=kin11_k(iom)+sig
2235 51680 : kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie)
2236 51680 : kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie)
2237 : kin22_k(iom)=kin22_k(iom) + &
2238 51680 : & sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie)
2239 : end if
2240 435600 : do l1=1,3
2241 1414400 : do l2=1,3
2242 : cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(l1,l2,iband,jband)&
2243 1305600 : & *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)
2244 : end do
2245 : end do
2246 :
2247 : end do
2248 :
2249 : ! Sumrule start
2250 420 : if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then
2251 : np_sum_k1=np_sum_k1 -dhdk2_g(iband,jband)&
2252 380 : & *(diff_occ)/(eig0_k(iband)-eig0_k(jband))
2253 : else
2254 20 : np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(iband,jband)
2255 : end if
2256 :
2257 :
2258 : ! end loop over band
2259 : end do
2260 21 : socc_k=socc_k+occ_k(iband)
2261 : end do
2262 :
2263 273 : do iom=1,mom
2264 272 : kin11(iom)=kin11(iom)+wtk(ikpt)*kin11_k(iom)
2265 272 : kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
2266 272 : kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
2267 272 : kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
2268 1089 : do l1=1,3
2269 3536 : do l2=1,3
2270 3264 : cond_kg(iom,l1,l2)=cond_kg(iom,l1,l2)+wtk(ikpt)*cond_nd(iom,l1,l2)
2271 : end do
2272 : end do
2273 : end do
2274 :
2275 1 : np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2)
2276 1 : socc=socc+wtk(ikpt)*socc_k
2277 :
2278 : ! Validity limit
2279 1 : deltae=deltae+(eig0_k(nband_k)-fermie)
2280 :
2281 1 : bd2tot_index=bd2tot_index+2*nband_k**2
2282 1 : bdtot_index=bdtot_index+nband_k
2283 1 : ABI_FREE(eig0_k)
2284 1 : ABI_FREE(eig1_k)
2285 1 : ABI_FREE(occ_k)
2286 1 : ABI_FREE(doccde_k)
2287 1 : ABI_FREE(dhdk2_r)
2288 2 : ABI_FREE(dhdk2_g)
2289 : ! End loop over k
2290 : end do
2291 :
2292 1 : write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three,socc
2293 : write(std_out,'(a,f10.5,a,f10.5,a)')&
2294 2 : & ' Emax-Efermi =',deltae/dble(nkpt),' Ha',deltae/dble(nkpt)*Ha_eV,' eV'
2295 :
2296 : !End loop over spins
2297 : end do
2298 :
2299 2461 : cond_kg=cond_kg*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
2300 :
2301 :
2302 : !Check that new output file does NOT exist
2303 : !Keep this line : prevent silly (compiler ?) bug on HP 8000
2304 1 : write(std_out,*)' conducti : call isfile '
2305 : !
2306 1 : if (open_file(trim(filnam_out)//'_tens',msg,newunit=tens_unt,form='formatted',action="write")/=0) then
2307 0 : ABI_ERROR(msg)
2308 : end if
2309 1 : if (open_file(trim(filnam_out)//'_Lij',msg,newunit=lij_unt,form='formatted',action="write")/=0) then
2310 0 : ABI_ERROR(msg)
2311 : end if
2312 1 : write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22'
2313 :
2314 1 : if (open_file(trim(filnam_out)//'_sig',msg,newunit=sig_unt,form='formatted',action="write")/=0) then
2315 0 : ABI_ERROR(msg)
2316 : end if
2317 1 : write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1'
2318 :
2319 1 : if (open_file(trim(filnam_out)//'_Kth',msg,newunit=kth_unt,form='formatted',action="write")/=0) then
2320 0 : ABI_ERROR(msg)
2321 : end if
2322 : write(kth_unt,'(a)')&
2323 1 : & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)'
2324 :
2325 1 : if (open_file(trim(filnam_out)//'.out',msg,newunit=ocond_unt,form='formatted',action="write")/=0) then
2326 0 : ABI_ERROR(msg)
2327 : end if
2328 1 : write(ocond_unt,'(a)' )' Conducti output file:'
2329 1 : write(ocond_unt,'(a)' )' Contains all results produced by conducti utility'
2330 1 : write(ocond_unt,'(a)' )' '
2331 1 : write(ocond_unt,'(a)')' # omega(ua) cond(ua) thermal cond(ua) thermopower(ua)'
2332 :
2333 : !Keep this line : prevent silly (compiler ?) bug on HP 8000
2334 1 : write(std_out,*)' conducti : after call isfile '
2335 :
2336 : !Compute thermal conductivity and thermopower
2337 273 : do iom=1,mom
2338 272 : oml=oml1(iom)
2339 272 : kin11(iom)=kin11(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
2340 272 : kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
2341 272 : kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
2342 272 : kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
2343 272 : if (dabs(kin11(iom))<10.0d-20) kin11(iom)=zero
2344 272 : Kth(iom)=kin22(iom)
2345 272 : Stp(iom)=zero
2346 272 : if(kin11(iom)/=zero) then
2347 107 : Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/kin11(iom))
2348 107 : Stp(iom)=kin12(iom)/(kin11(iom)*Tatm)
2349 : end if
2350 272 : if (dabs(Kth(iom))<10.0d-20) Kth(iom)=zero
2351 272 : if (dabs(Stp(iom))<10.0d-20) Stp(iom)=zero
2352 272 : if (dabs(kin12(iom))<10.0d-20) kin12(iom)=zero
2353 272 : if (dabs(kin21(iom))<10.0d-20) kin21(iom)=zero
2354 272 : if (dabs(kin22(iom))<10.0d-20) kin22(iom)=zero
2355 :
2356 272 : write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9
2357 272 : write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,kin11(iom),kin11(iom)*Ohmcm
2358 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
2359 273 : write(ocond_unt,'(1f12.5,3es22.12)') oml,kin11(iom),Kth(iom),Stp(iom)
2360 : end do
2361 :
2362 1 : write(tens_unt,'(a)' )' Conductivity file '
2363 1 : write(tens_unt,'(a)' )' ----------------- '
2364 1 : write(tens_unt,'(a)' )' Contain first the full conductivity tensor, for the desired set of energies,'
2365 1 : write(tens_unt,'(a)' )' then, the three principal values, for the desired set of energies'
2366 1 : write(tens_unt,'(a)' )' (note that eigenvalues are not directly associated with xx,yy,zz)'
2367 1 : write(tens_unt,'(a)' )' '
2368 :
2369 1 : write(ocond_unt,'(a)' )' '
2370 1 : write(ocond_unt,'(a)' )' full conductivity tensor, for the desired set of energies'
2371 1 : write(ocond_unt,'(a)' )' then, the three principal values, for the desired set of energies:'
2372 :
2373 273 : do iom=1,mom
2374 272 : oml=oml1(iom)*Ha_eV
2375 272 : write(tens_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
2376 272 : write(ocond_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
2377 1089 : do l1=1,3
2378 3264 : write(tens_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
2379 3536 : write(ocond_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
2380 : end do
2381 : end do
2382 :
2383 : !Diagonalizing the conductivity matrix for sigma_xx,sigma_yy,sigma_zz
2384 273 : cond_kg_xx=0d0
2385 273 : cond_kg_yy=0d0
2386 273 : cond_kg_zz=0d0
2387 : !trace=0d0 ! Used for checking with the original version of the code
2388 273 : do iom=1,mom
2389 272 : oml=oml1(iom)*Ha_eV
2390 272 : cond_kg_w=0d0
2391 1088 : do l1=1,3
2392 3536 : do l2=1,3
2393 3264 : cond_kg_w(l1,l2)=cond_kg(iom,l1,l2)
2394 : end do
2395 : end do
2396 272 : call jacobi(cond_kg_w,3,3,eig_cond,z,nrot)
2397 :
2398 : ! When the value is too small, set it to zero before printing
2399 272 : if(abs(eig_cond(1))<tol10)eig_cond(1)=zero
2400 272 : if(abs(eig_cond(2))<tol10)eig_cond(2)=zero
2401 272 : if(abs(eig_cond(3))<tol10)eig_cond(3)=zero
2402 :
2403 272 : cond_kg_xx(iom)=eig_cond(1)
2404 272 : cond_kg_yy(iom)=eig_cond(2)
2405 273 : cond_kg_zz(iom)=eig_cond(3)
2406 : ! trace(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
2407 : end do
2408 :
2409 : !DEBUG Keep this line : prevent silly (compiler ?) bug on HP 8000
2410 : !write(std_out,*)' conducti : after open '
2411 : !ENDDEBUG
2412 :
2413 1 : write(tens_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
2414 1 : write(tens_unt,'(a)')' '
2415 1 : write(tens_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)'
2416 :
2417 1 : write(ocond_unt,'(a)')' '
2418 1 : write(ocond_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
2419 1 : write(ocond_unt,'(a)')' '
2420 1 : write(ocond_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)'
2421 :
2422 :
2423 273 : do iom=1,mom
2424 272 : cond_tot(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
2425 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)
2426 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)
2427 : end do
2428 :
2429 1 : write(tens_unt,*)
2430 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'
2431 1 : write(ocond_unt,*)
2432 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'
2433 :
2434 273 : do iom=1,mom
2435 272 : oml=oml1(iom)*Ha_eV
2436 272 : cond_tot(iom)=cond_tot(iom)*Ohmcm
2437 272 : cond_kg_xx(iom)=cond_kg_xx(iom)*Ohmcm
2438 272 : cond_kg_yy(iom)=cond_kg_yy(iom)*Ohmcm
2439 272 : cond_kg_zz(iom)=cond_kg_zz(iom)*Ohmcm
2440 272 : write(tens_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
2441 273 : write(ocond_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
2442 : end do
2443 :
2444 : !Calculate the imaginary part of the conductivity (principal value)
2445 : !+derived optical properties.
2446 1 : call msig(kin11,mom,oml1,filnam_out,zero,0)
2447 :
2448 1 : close(tens_unt)
2449 1 : close(lij_unt)
2450 1 : close(sig_unt)
2451 1 : close(kth_unt)
2452 1 : close(ocond_unt)
2453 :
2454 1 : ABI_FREE(nband)
2455 1 : ABI_FREE(oml1)
2456 1 : ABI_FREE(occ)
2457 1 : ABI_FREE(eigen11)
2458 1 : ABI_FREE(eigen12)
2459 1 : ABI_FREE(eigen13)
2460 1 : ABI_FREE(eigen0)
2461 1 : ABI_FREE(doccde)
2462 1 : ABI_FREE(wtk)
2463 1 : ABI_FREE(cond_nd)
2464 1 : ABI_FREE(cond_kg)
2465 1 : ABI_FREE(cond_kg_cart)
2466 1 : ABI_FREE(cond_kg_xx)
2467 1 : ABI_FREE(cond_kg_yy)
2468 1 : ABI_FREE(trace)
2469 1 : ABI_FREE(cond_kg_zz)
2470 1 : ABI_FREE(cond_tot)
2471 1 : ABI_FREE(kin11)
2472 1 : ABI_FREE(kin22)
2473 1 : ABI_FREE(kin12)
2474 1 : ABI_FREE(kin21)
2475 1 : ABI_FREE(kin11_k)
2476 1 : ABI_FREE(kin22_k)
2477 1 : ABI_FREE(kin12_k)
2478 1 : ABI_FREE(kin21_k)
2479 1 : ABI_FREE(Stp)
2480 1 : ABI_FREE(Kth)
2481 1 : call hdr%free()
2482 :
2483 5 : end subroutine conducti_nc
2484 : !!***
2485 :
2486 : !----------------------------------------------------------------------
2487 :
2488 : !!****f* m_conducti/msig
2489 : !! NAME
2490 : !! msig
2491 : !!
2492 : !! FUNCTION
2493 : !! This program computes the elements of the optical frequency dependent
2494 : !! conductivity tensor and the conductivity along the three principal axes
2495 : !! from the Kubo-Greenwood formula for PAW formalism
2496 : !!
2497 : !! INPUTS
2498 : !! fcti(npti)= conductivity, as calculated in conducti
2499 : !! npti= number of points to calculate conductivity
2500 : !! xi(npti)= energies where the conductivity is calculated
2501 : !!
2502 : !! OUTPUT
2503 : !! no output, only files
2504 : !!
2505 : !! NOTES
2506 : !! this program calculates the imaginary part of the conductivity (principal value)
2507 : !! +derived optical properties.
2508 : !!
2509 : !! SOURCE
2510 :
2511 3 : subroutine msig(fcti,npti,xi,filnam_out_sig,phi,au_units)
2512 :
2513 : !Arguments -----------------------------------
2514 : !scalars
2515 : integer,intent(in) :: npti, au_units
2516 : !arrays
2517 : real(dp),intent(in) :: fcti(npti),xi(npti)
2518 : character(len=fnlen),intent(in) :: filnam_out_sig
2519 : real(dp),intent(in) :: phi
2520 :
2521 : !Local variables-------------------------------
2522 : !scalars
2523 : integer :: ii,ip,eps_unt,abs_unt
2524 : real(dp),parameter :: del=0.001_dp
2525 : real(dp) :: dx,eps1,eps2,komega,pole,refl_s,refl_p,sigma2,xsum,ff,ffp,ffpp,abso,sigma1
2526 : character(len=500) :: msg
2527 : !arrays
2528 3 : real(dp),allocatable :: fct(:)
2529 : real(dp) :: xx1,xx2,xx,nomega
2530 : complex(dp) :: epsc,cos_phi,sin_phi,sqroot,crefl_s,crefl_p
2531 : ! *********************************************************************************
2532 :
2533 :
2534 3 : write(std_out,'(2a)')ch10,'Calculate the principal value and related optical properties'
2535 3 : write(std_out,'(a)')'Use default value for delta interval: del=1e-3'
2536 :
2537 3 : if (open_file(trim(filnam_out_sig)//'_eps',msg,newunit=eps_unt,status='replace',action="write")/=0) then
2538 0 : ABI_ERROR(msg)
2539 : end if
2540 3 : if(au_units==0) then
2541 3 : write(eps_unt,'(a)')'#energy (eV),sigma_1(Ohm-1cm-1),sigma_2(Ohm-1cm-1),epsilon_1(cgs),epsilon_2(cgs)'
2542 : else
2543 0 : write(eps_unt,'(a)')'#energy (Ha),sigma_1(au),sigma_2(au),epsilon_1(au),epsilon_2(au)'
2544 : endif
2545 :
2546 3 : if (open_file(trim(filnam_out_sig)//'_abs',msg,newunit=abs_unt,status='replace',action="write")/=0) then
2547 0 : ABI_ERROR(msg)
2548 : end if
2549 3 : if(au_units==0) then
2550 3 : write(abs_unt,'(a)')'#energy(eV),nomega,komega,refl. s,refl. p,abso.(cm-1)'
2551 : else
2552 0 : write(abs_unt,'(a)')'#energy(Ha),nomega,komega,refl. s,refl. p,abso.(au)'
2553 : endif
2554 :
2555 9 : ABI_MALLOC(fct,(npti))
2556 :
2557 : !loop on the initial energy grid
2558 2275 : do ip=1,npti
2559 : !!! Taylor expansion up to second order of fcti at xx=xi(ii) for each ii
2560 : !!! Then, in each discretization interval, the integral can be performed analytically
2561 2272 : pole=xi(ip)
2562 2272 : dx=(xi(npti)-xi(1))/dble(npti-1)
2563 2272 : xsum=zero
2564 2076256 : do ii=1,npti
2565 2073984 : xx=xi(ii)
2566 2076256 : fct(ii)=fcti(ii)*pole/(xx+pole)
2567 : enddo
2568 2076256 : do ii=1,npti
2569 2073984 : xx=xi(ii)
2570 2073984 : ff=fct(ii)
2571 2073984 : ffp=zero
2572 2073984 : ffpp=zero
2573 2073984 : if(ii<=npti-3) then
2574 2067168 : ffp=(four*fct(ii+1)-three*fct(ii)-fct(ii+2))/dx/two
2575 : ffpp=(-fct(ii+3)+four*fct(ii+2)-five*fct(ii+1)+&
2576 2067168 : & two*fct(ii))/dx/dx
2577 : endif
2578 2073984 : xx2=xx+half*dx
2579 2073984 : xx1=zero
2580 2073984 : if(ii>1) xx1=xx-half*dx
2581 : xsum=xsum+ff*log(abs((xx2-pole)/(xx1-pole)))+ffp*(xx2-xx1+(pole-xx)*log(abs((xx2-pole)/(xx1-pole))))+&
2582 2076256 : & half*ffpp*((xx-pole)**2*log(abs((xx2-pole)/(xx1-pole)))+(xx2**2-xx1**2+(two*pole-four*xx)*(xx2-xx1))/two)
2583 : enddo
2584 2272 : if(pole<tol3) xsum=zero
2585 :
2586 : ! Calculate the derivated optical quantities and output the value
2587 2272 : sigma2=(-two/pi)*xsum
2588 2272 : eps1=one-(four_pi*sigma2/(pole))
2589 2272 : eps2=four*fcti(ip)*pi/(pole)
2590 :
2591 : ! A special treatment of the case where eps2 is very small compared to eps1 is needed
2592 2272 : if(eps2**2 > eps1**2 * tol12)then
2593 1172 : nomega=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
2594 1172 : komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
2595 1172 : abso=four_pi*fcti(ip)/(nomega*Speed_Light)
2596 1100 : else if(eps1>zero)then
2597 921 : nomega=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
2598 921 : komega=half*abs(eps2/sqrt(eps1))
2599 921 : abso=four_pi*fcti(ip)/(nomega*Speed_Light)
2600 179 : else if(eps1<zero)then
2601 179 : nomega=half*abs(eps2/sqrt(-eps1))
2602 179 : komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
2603 179 : abso=two*sqrt(-eps1)*pole/(Speed_Light)
2604 : end if
2605 :
2606 2272 : epsc=cmplx(eps1,eps2,kind=dp)
2607 2272 : cos_phi=cmplx(cos(phi),kind=dp)
2608 2272 : sin_phi=cmplx(sin(phi),kind=dp)
2609 2272 : sqroot=sqrt(epsc-sin_phi*sin_phi)
2610 2272 : crefl_s=(cos_phi-sqroot)/(cos_phi+sqroot)
2611 2272 : crefl_p=(epsc*cos_phi-sqroot)/(epsc*cos_phi+sqroot)
2612 :
2613 2272 : refl_s=real(crefl_s*conjg(crefl_s))
2614 2272 : refl_p=real(crefl_p*conjg(crefl_p))
2615 :
2616 2272 : sigma1=fcti(ip)
2617 2272 : if(au_units==0) then
2618 2272 : pole=pole*Ha_eV
2619 2272 : sigma1=sigma1*Ohmcm
2620 2272 : sigma2=sigma2*Ohmcm
2621 2272 : abso=abso*Ohmcm*Speed_Light/(Speed_Light_SI*four_pi*eps0)
2622 : endif
2623 2272 : write(eps_unt,'(5e18.10)') pole,sigma1,sigma2,eps1,eps2
2624 2275 : write(abs_unt,'(6e18.10)') pole,nomega,komega,refl_s,refl_p,abso
2625 :
2626 : end do
2627 :
2628 3 : close(eps_unt)
2629 3 : close(abs_unt)
2630 :
2631 3 : ABI_FREE(fct)
2632 :
2633 3 : end subroutine msig
2634 : !!***
2635 :
2636 : end module m_conducti
2637 : !!***
|
