Branch data Line data Source code
1 : : !!****m* ABINIT/m_conducti
2 : : !! NAME
3 : : !! m_conducti
4 : : !!
5 : : !! FUNCTION
6 : : !! This program computes the elements of the optical frequency dependent
7 : : !! conductivity tensor and the conductivity along the three principal axes
8 : : !! from the Kubo-Greenwood formula.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2002-2026 ABINIT group (VRecoules, PGhosh, SMazevet, SM, SVinko, NBrouwer)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : 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 : : !!***
|
