Line data Source code
1 : !!****m* ABINIT/m_thmeig
2 : !! NAME
3 : !! m_thmeig
4 : !!
5 : !! FUNCTION
6 : !! Calculate thermal corrections to the eigenvalues.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (PB, XG, GA)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_thmeig
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_tetrahedron
27 : use m_htetra
28 : use m_errors
29 : use m_ddb
30 : use m_ddb_hdr
31 : use m_xmpi
32 : use m_sort
33 :
34 : use m_geometry, only : mkrdim, xred2xcart, metric
35 : use m_symfind, only : symfind, symlatt
36 : use m_matrix, only : mati3inv, matr3inv
37 : use m_symtk, only : symatm
38 : use m_crystal, only : crystal_t
39 : use m_io_tools, only : open_file
40 : use m_dynmat, only : asria_corr, dfpt_phfrq
41 : use m_anaddb_dataset, only : anaddb_dataset_type
42 : use m_pawtab, only : pawtab_type,pawtab_nullify,pawtab_free
43 : use m_kpts, only : getkgrid
44 :
45 : implicit none
46 :
47 : private
48 : !!***
49 :
50 : public :: thmeig
51 : !!***
52 :
53 : contains
54 : !!***
55 :
56 : !!****f* m_thmeig/thmeig
57 : !! NAME
58 : !! thmeig
59 : !!
60 : !! FUNCTION
61 : !! This routine calculates the thermal corrections to the eigenvalues.
62 : !! The output is this quantity for the input k point.
63 : !!
64 : !! INPUTS
65 : !! elph_base_name = root filename for outputs
66 : !! eig2_filnam = name of the eig2 database file
67 : !! comm=MPI communicator
68 : !!
69 : !! OUTPUT
70 : !!
71 : !! SOURCE
72 :
73 7 : subroutine thmeig(inp, ddb, crystal, iout, natom, mpert, msize, d2asr, comm)
74 :
75 : !Arguments ------------------------------------
76 : !scalars
77 : integer,intent(inout) :: natom
78 : integer,intent(in) :: mpert,msize
79 : integer,intent(in) :: comm
80 : integer,intent(in) :: iout
81 : type(crystal_t), intent(inout) :: crystal
82 : type(anaddb_dataset_type),intent(inout) :: inp
83 : type(ddb_type),intent(inout) :: ddb
84 : !arrays
85 : real(dp),intent(inout) :: d2asr(2,3,natom,3,natom)
86 :
87 :
88 : !Local variables-------------------------------
89 : !scalars
90 : integer,parameter :: msppol=2,master=0,bcorr0=0
91 : integer :: msym
92 : integer :: nkpt,mband,ntypat
93 : integer :: usepaw,natifc
94 : integer :: nsym,occopt,nblok2
95 : integer :: ntemper,telphint,thmflag
96 : integer :: brav,chksymbreak,found,gqpt,iatom1,iatom2,iband,iblok,iblok2,idir1,idir2,ii,ikpt,ilatt,imod,index
97 : integer :: iomega,iqpt,iqpt1,iqpt2,iqpt2_previous,iqpt3,iscf_fake,itemper
98 : integer :: mpert_eig2,msize2,nene,ng2f,nqshft,nsym_new,unit_g2f,nqpt,nqpt_computed,qptopt,rftyp
99 : !integer :: mqpt,nqpt2,option
100 : integer :: unit_phdos,unitout
101 : integer :: isym
102 : integer :: nptsym,use_inversion
103 : integer :: ierr
104 : real(dp) :: ucvol
105 : real(dp) :: g2fsmear,temperinc,tempermin
106 : real(dp) :: bosein,deltaene,det,domega,enemax,enemin,fact2i,fact2r,factr
107 : real(dp) :: gaussfactor,gaussprefactor,gaussval,invdet,omega,omega_max,omega_min,qnrm,qptrlen
108 : real(dp) :: rcvol,tmp,tol,vec1i,vec1r,vec2i,vec2r,veci,vecr,xx
109 : real(dp) :: tolsym,tolsym8 !new
110 : character(len=500) :: message
111 : character(len=fnlen) :: outfile, elph_base_name, eig2_filnam
112 7 : type(ddb_type) :: ddb_eig2
113 364 : type(ddb_hdr_type) :: ddb_hdr
114 : !arrays
115 : ! FIXME now these must be allocated
116 : integer :: ngqpt(9),qptrlatt(3,3),rfelfd(4),rfphon(4),rfstrs(4),vacuum(3)
117 : integer :: bravais(11)
118 7 : integer,allocatable :: typat(:)
119 7 : integer,allocatable :: symrel(:,:,:),symrec(:,:,:)
120 7 : integer,allocatable :: indsym(:,:,:)
121 7 : integer,allocatable :: indqpt(:)
122 7 : integer,allocatable :: symafm(:),symafm_new(:)
123 7 : integer,allocatable :: carflg_eig2(:,:,:,:)
124 7 : integer,allocatable :: ptsymrel(:,:,:),symrel_new(:,:,:)
125 : !integer,allocatable :: symrec_new(:,:,:)
126 : real(dp) :: rprim(3,3),gprim(3,3),rmet(3,3),gmet(3,3)
127 : real(dp) :: acell(3)
128 : real(dp) :: diff_qpt(3)
129 : real(dp) :: gprimd(3,3),mesh(3,3)
130 : real(dp) :: qlatt(3,3),qphnrm(3),qpt_search(3,3)
131 : real(dp) :: rprimd(3,3),shiftq(3,MAX_NSHIFTK),tempqlatt(3)
132 : real(dp) :: dummy(0),dummy2(0,0)
133 7 : real(dp),allocatable :: xcart(:,:),xred(:,:)
134 14 : real(dp),allocatable :: amu(:),zion(:)
135 7 : real(dp),allocatable :: tnons(:,:)
136 7 : real(dp),allocatable :: deigi(:,:), deigr(:,:), multi(:,:), multr(:,:)
137 7 : real(dp),allocatable :: dwtermi(:,:), dwtermr(:,:)
138 7 : real(dp),allocatable :: slope(:,:,:),thmeigen(:,:,:),zeropoint(:,:,:)
139 7 : real(dp),allocatable :: displ(:)
140 7 : real(dp),allocatable :: dos_phon(:),dtweightde(:,:),d2cart(:,:)
141 7 : real(dp),allocatable :: eigvec(:,:,:,:),eigval(:,:),g2f(:,:,:),intweight(:,:,:)
142 7 : real(dp),allocatable :: indtweightde(:,:,:),tmpg2f(:,:,:),tmpphondos(:),total_dos(:),tweight(:,:)
143 7 : real(dp),allocatable :: phfreq(:,:)
144 7 : real(dp),allocatable :: eig2dGamma(:,:,:,:),kpnt(:,:,:)
145 7 : real(dp),allocatable :: dedni(:,:,:,:),dednr(:,:,:,:)
146 7 : real(dp),allocatable :: eigen_in(:)
147 7 : real(dp),allocatable :: qpt_full(:,:),qptnrm(:)
148 7 : real(dp),allocatable :: spqpt(:,:),tnons_new(:,:),spinat(:,:)
149 7 : real(dp),allocatable :: wghtq(:)
150 :
151 7 : type(t_tetrahedron) :: tetrahedra
152 : !type(htetra_t) :: tetrahedra
153 : character(len=80) :: errstr
154 :
155 : ! *********************************************************************
156 :
157 : ! Only master works for the time being
158 7 : if (xmpi_comm_rank(comm) /= master) return
159 :
160 574 : write(message,'(83a)') ch10,('=',ii=1,80),ch10,&
161 581 : & ' Computation of the electron-phonon changes to the electronic eigenenergies '
162 7 : call wrtout(ab_out,message,'COLL')
163 7 : call wrtout(std_out,message,'COLL')
164 :
165 :
166 : !=========================================================================
167 : !0) Initializations
168 : !=========================================================================
169 :
170 : !GA: TODO Perhaps those should be created at initialization of inp
171 7 : elph_base_name = trim(inp%prefix_outdata)//"_ep"
172 7 : eig2_filnam = inp%filename_eigr2d
173 :
174 :
175 7 : g2fsmear = inp%a2fsmear
176 :
177 7 : telphint = inp%telphint
178 7 : temperinc = inp%temperinc
179 7 : tempermin = inp%tempermin
180 7 : thmflag = inp%thmflag
181 :
182 7 : ntemper = inp%ntemper
183 7 : natifc = inp%natifc
184 :
185 : !Open Derivative DataBase then r/w Derivative DataBase preliminary information.
186 :
187 7 : write(std_out, '(a)' ) '- thmeig: Initialize the second-order electron-phonon file with name :'
188 7 : write(std_out, '(a,a)' )'- ',trim(eig2_filnam)
189 :
190 7 : call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self)
191 :
192 7 : mband = ddb_hdr%mband * ddb_hdr%nsppol
193 7 : nkpt = ddb_hdr%nkpt
194 7 : ntypat = ddb_hdr%ntypat
195 :
196 7 : msym = ddb_hdr%msym
197 7 : nblok2 = ddb_hdr%nblok
198 7 : nsym = ddb_hdr%nsym
199 7 : occopt = ddb%occopt
200 7 : usepaw = ddb_hdr%usepaw
201 :
202 21 : ABI_MALLOC(typat, (natom))
203 21 : ABI_MALLOC(zion, (ntypat))
204 14 : ABI_MALLOC(amu, (ntypat))
205 :
206 21 : ABI_MALLOC(xcart,(3,natom))
207 14 : ABI_MALLOC(xred,(3,natom))
208 :
209 21 : ABI_MALLOC(symafm, (nsym))
210 14 : ABI_MALLOC(spinat,(3,natom))
211 :
212 21 : ABI_MALLOC(symrel, (3,3,nsym))
213 14 : ABI_MALLOC(symrec, (3,3,nsym))
214 21 : ABI_MALLOC(tnons, (3,nsym))
215 28 : ABI_MALLOC(indsym, (4,nsym,natom))
216 :
217 28 : ABI_MALLOC(deigi, (mband,nkpt))
218 21 : ABI_MALLOC(deigr, (mband,nkpt))
219 21 : ABI_MALLOC(dwtermi, (mband,nkpt))
220 21 : ABI_MALLOC(dwtermr, (mband,nkpt))
221 21 : ABI_MALLOC(multi, (mband,nkpt))
222 21 : ABI_MALLOC(multr, (mband,nkpt))
223 28 : ABI_MALLOC(slope, (2,mband,nkpt))
224 21 : ABI_MALLOC(thmeigen, (2,mband,nkpt))
225 21 : ABI_MALLOC(zeropoint, (2,mband,nkpt))
226 :
227 : !At present, only atom-type perturbations are allowed for eig2 type matrix elements.
228 7 : mpert_eig2=natom
229 7 : msize2=3*mpert_eig2*3*mpert_eig2
230 :
231 7 : ddb_eig2%nsppol = ddb_hdr%nsppol
232 7 : call ddb_eig2%malloc(msize2,nblok2,natom,ntypat,mpert_eig2,nkpt,mband)
233 :
234 35 : ABI_MALLOC(eig2dGamma,(2,msize2,mband,nkpt))
235 :
236 28 : ABI_MALLOC(eigvec,(2,3,natom,3*natom))
237 28 : ABI_MALLOC(phfreq,(3*natom,ddb%nblok))
238 :
239 : !amu = ddb%amu
240 14 : amu(:) = ddb_hdr%amu(1:ntypat)
241 26 : typat(:) = ddb_hdr%typat(1:natom)
242 14 : zion(:) = ddb_hdr%zion(1:ntypat)
243 98 : symrel(:,:,1:nsym) = ddb_hdr%symrel(:,:,1:nsym)
244 35 : tnons(:,1:nsym) = ddb_hdr%tnons(:,1:nsym)
245 :
246 83 : xred(:,:) = ddb_hdr%xred(:,:)
247 :
248 14 : symafm(:) = ddb_hdr%symafm(1:nsym)
249 83 : spinat(:,:) = ddb_hdr%spinat(:,1:natom)
250 :
251 : !symrel = ddb_hdr%symrel ! out
252 : !tnons = ddb_hdr%tnons ! out
253 :
254 : !acell = ddb%acell
255 : !natom = ddb_hdr%natom
256 28 : acell = ddb_hdr%acell
257 91 : rprim = ddb_hdr%rprim
258 :
259 : !Compute different matrices in real and reciprocal space, also
260 : !checks whether ucvol is positive.
261 7 : call mkrdim(acell,rprim,rprimd)
262 7 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
263 :
264 : !Obtain reciprocal space primitive transl g from inverse trans of r
265 7 : call matr3inv(rprim,gprim)
266 :
267 : !Generate atom positions in cartesian coordinates
268 7 : call xred2xcart(natom,rprimd,xcart,xred)
269 :
270 : !Transposed inversion of the symmetry matrices, for use in
271 : !the reciprocal space
272 14 : do isym=1,nsym
273 14 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
274 : end do
275 :
276 : !SYMATM generates for all the atoms and all the symmetries, the atom
277 : !on which the referenced one is sent and also the translation bringing
278 : !back this atom to the referenced unit cell
279 7 : tolsym8=tol8
280 7 : call symatm(indsym,natom,nsym,symrec(:,:,1:nsym),tnons(:,1:nsym),tolsym8,typat,xred)
281 :
282 91391 : eig2dGamma(:,:,:,:)=zero
283 :
284 35 : ABI_MALLOC(carflg_eig2,(3,mpert_eig2,3,mpert_eig2))
285 21 : ABI_MALLOC(kpnt,(3,nkpt,1))
286 :
287 : ! Copy a bunch of stuff back into crystal (to retain old behavior)
288 : ! TODO comment these: doesnt make a difference
289 90 : crystal%xcart = xcart
290 7 : crystal%ucvol = ucvol
291 21 : crystal%zion = zion
292 91 : crystal%gmet = gmet
293 91 : crystal%rmet = rmet
294 7 : crystal%nsym = nsym
295 105 : crystal%symrel = symrel
296 105 : crystal%symrec = symrec
297 42 : crystal%tnons = tnons
298 128 : crystal%indsym = indsym
299 :
300 :
301 : !=========================================================================
302 : !1) Take care of the Gamma point for thmflag=3, 5 or 7
303 : !=========================================================================
304 :
305 7 : if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
306 7 : found=0
307 7 : do iblok2=1,nblok2
308 :
309 7 : call ddb_eig2%read_d2eig(ddb_hdr, iblok2, iblok2)
310 :
311 : qnrm = ddb_eig2%qpt(1,iblok2)*ddb_eig2%qpt(1,iblok2)+ &
312 : & ddb_eig2%qpt(2,iblok2)*ddb_eig2%qpt(2,iblok2)+ &
313 7 : & ddb_eig2%qpt(3,iblok2)*ddb_eig2%qpt(3,iblok2)
314 7 : if(qnrm < DDB_QTOL) then
315 91391 : eig2dGamma(:,:,:,:) = ddb_eig2%eig2dval(:,:,:,:)
316 7 : gqpt=iblok2
317 7 : write(std_out,*)'-thmeig: found Gamma point in EIG2 DDB, blok number ',iblok2
318 : found=1
319 : exit
320 : end if
321 : end do
322 :
323 : if(found==0)then
324 : write(message,'(a,i3,2a)')&
325 0 : & 'Was unable to find the blok for Gamma point in EIG2 DDB file, while thmflag= ',thmflag,ch10,&
326 0 : & 'Action: compute the contribution from Gamma, and merge it in your EIG2 DDB file.'
327 0 : ABI_ERROR(message)
328 : end if
329 :
330 : ! Put eig2dGamma in cartesian coordinates
331 : call carttransf(ddb_eig2%flg,eig2dGamma,carflg_eig2,gprimd,gqpt,mband,&
332 7 : & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd)
333 :
334 : end if
335 :
336 7 : call ddb_hdr%close()
337 :
338 : !=========================================================================
339 : !2) Calculation of dE(n,k)/dn(Q,j) : consider all q and modes
340 : !=========================================================================
341 :
342 7 : if(thmflag==3 .or. thmflag==4)then
343 :
344 :
345 : ! Use the first list of q wavevectors
346 4 : nqpt=inp%nph1l
347 12 : ABI_MALLOC(spqpt,(3,nqpt))
348 12 : do iqpt=1,inp%nph1l
349 36 : spqpt(:,iqpt)=inp%qph1l(:,iqpt)/inp%qnrml1(iqpt)
350 : end do
351 12 : ABI_MALLOC(wghtq,(nqpt))
352 12 : wghtq(:)=one/nqpt
353 :
354 3 : else if(thmflag>=5 .and. thmflag<=8)then
355 :
356 : ! Generates the q point grid
357 12 : ngqpt(1:3)=inp%ngqpt(1:3)
358 3 : nqshft=inp%nqshft
359 3 : qptrlatt(:,:)=0
360 3 : qptrlatt(1,1)=ngqpt(1)
361 3 : qptrlatt(2,2)=ngqpt(2)
362 3 : qptrlatt(3,3)=ngqpt(3)
363 :
364 9 : ABI_MALLOC(ptsymrel,(3,3,msym))
365 9 : ABI_MALLOC(symafm_new,(msym))
366 6 : ABI_MALLOC(symrel_new,(3,3,msym))
367 9 : ABI_MALLOC(tnons_new,(3,msym))
368 3 : if(thmflag==7 .or. thmflag==8) then
369 : ! Re-generate symmetry operations from the lattice and atomic coordinates
370 2 : tolsym=tol8
371 2 : call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym)
372 2 : use_inversion=1
373 : call symfind(gprimd,msym,natom,nptsym,1,nsym_new,0,&
374 2 : & ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred)
375 2 : write(std_out,*)' thmeig : found ',nsym_new,' symmetries ',ch10
376 2 : qptopt=1
377 : else
378 1 : nsym_new=1
379 13 : symrel_new(:,:,1)=0 ; symrel_new(1,1,1)=1 ; symrel_new(2,2,1)=1 ; symrel_new(3,3,1)=1
380 4 : tnons_new(:,1)=zero
381 1 : symafm_new(1)=1
382 1 : qptopt=3
383 : end if
384 :
385 : ! GA: This is useless
386 3 : brav=inp%brav
387 3 : if(abs(brav)/=1)then
388 0 : message = ' The possibility to have abs(brav)/=1 for thmeig was disabled.'
389 0 : ABI_ERROR(message)
390 : end if
391 :
392 : ! Prepare to compute the q-point grid in the ZB or IZB
393 3 : iscf_fake=5 ! Need the weights
394 3 : chksymbreak=0
395 3 : vacuum=0
396 27 : shiftq(:,1:nqshft)=inp%q1shft(:,1:nqshft)
397 : ! Compute the final number of q points
398 : call getkgrid(chksymbreak,0,iscf_fake,dummy2,qptopt,qptrlatt,qptrlen,&
399 : & nsym_new,0,nqpt,nqshft,nsym_new,rprimd,&
400 3 : & shiftq,symafm_new,symrel_new,vacuum,dummy)
401 9 : ABI_MALLOC(spqpt,(3,nqpt))
402 9 : ABI_MALLOC(wghtq,(nqpt))
403 : call getkgrid(chksymbreak,iout,iscf_fake,spqpt,qptopt,qptrlatt,qptrlen,&
404 : & nsym_new,nqpt,nqpt_computed,nqshft,nsym_new,rprimd,&
405 3 : & shiftq,symafm_new,symrel_new,vacuum,wghtq)
406 :
407 3 : ABI_FREE(ptsymrel)
408 3 : ABI_FREE(symafm_new)
409 3 : ABI_FREE(symrel_new)
410 6 : ABI_FREE(tnons_new)
411 :
412 : end if
413 :
414 7 : call ddb_hdr%free()
415 :
416 :
417 7 : write(message,'(a,a)')ch10,' thmeig : list of q wavevectors, with integration weights '
418 7 : call wrtout(ab_out,message,'COLL')
419 7 : call wrtout(std_out,message,'COLL')
420 28 : do iqpt=1,nqpt
421 21 : write(message,'(i6,3es16.6,es20.6)')iqpt,spqpt(:,iqpt),wghtq(iqpt)
422 21 : call wrtout(ab_out,message,'COLL')
423 28 : call wrtout(std_out,message,'COLL')
424 : end do
425 :
426 : if(.not.allocated(indqpt))then
427 21 : ABI_MALLOC(indqpt,(nqpt))
428 : end if
429 42 : ABI_MALLOC(dedni,(mband,nkpt,3*natom,nqpt))
430 35 : ABI_MALLOC(dednr,(mband,nkpt,3*natom,nqpt))
431 21 : ABI_MALLOC(eigen_in,(nqpt))
432 21 : ABI_MALLOC(qpt_full,(3,nqpt))
433 14 : ABI_MALLOC(qptnrm,(nqpt))
434 :
435 7246 : dednr(:,:,:,:) = zero
436 7246 : dedni(:,:,:,:) = zero
437 :
438 : !!Prepare the reading of the EIG2 files
439 7 : call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self, msym=msym)
440 7 : call ddb_hdr%free()
441 :
442 : !iqpt2 will be the index of the q point bloks inside the EIG2 file
443 7 : iqpt2=0
444 :
445 : !Sum on all phonon wavevectors and modes
446 28 : do iqpt=1,nqpt
447 :
448 : ! Finding the target wavevector in DDB file
449 21 : qpt_search(:,:)=0.0d0
450 84 : qpt_search(:,1)=spqpt(:,iqpt)
451 84 : qphnrm(:)=one
452 63 : rfphon(1:2)=1
453 : ! NOTE : at present, no LO-TO splitting included !!!
454 63 : rfelfd(1:2)=0
455 63 : rfstrs(1:2)=0
456 21 : rftyp=1
457 :
458 21 : write(std_out,'(a,3es16.6)' )' Looking for spqpt=',qpt_search(:,1)
459 :
460 21 : call ddb%get_block(iblok,qpt_search,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
461 :
462 21 : if(iblok==0) then
463 : write(message,'(a,3es16.6,2a)')&
464 0 : & 'Was unable to find in DDB file, the blok for point ',spqpt(:,iqpt),ch10,&
465 0 : & 'Action: compute the contribution from this point, and merge it in your DDB file.'
466 0 : ABI_ERROR(message)
467 : end if
468 :
469 63 : ABI_MALLOC(d2cart,(2,msize))
470 : ! Copy the dynamical matrix in d2cart
471 111612 : d2cart(:,1:msize)=ddb%val(:,:,iblok)
472 :
473 : ! Eventually impose the acoustic sum rule based on previously calculated d2asr
474 : !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
475 21 : if (inp%asr==1 .or. inp%asr==2 .or. inp%asr==5) then
476 21 : call asria_corr(inp%asr,d2asr,d2cart,mpert,natom)
477 : end if
478 :
479 : ! Calculation of the eigenvectors and eigenvalues
480 : ! of the dynamical matrix
481 63 : ABI_MALLOC(displ,(2*3*natom*3*natom))
482 63 : ABI_MALLOC(eigval,(3,natom))
483 : call dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,&
484 : & mpert,msym,natom,nsym,ntypat,phfreq(:,iqpt),qphnrm(1),spqpt(:,iqpt),rprimd,inp%symdynmat,&
485 21 : & symrel,symafm,typat,ucvol)
486 21 : ABI_FREE(displ)
487 21 : ABI_FREE(eigval)
488 21 : ABI_FREE(d2cart)
489 :
490 :
491 : ! Read the next bloks to find the next q point.
492 21 : found=0 ; iqpt2_previous=iqpt2
493 21 : do while (iqpt2<nblok2)
494 21 : iqpt2=iqpt2+1
495 21 : call ddb_eig2%read_d2eig(ddb_hdr, iqpt2, iqpt2)
496 84 : diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt)
497 21 : if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then
498 : found=1
499 : exit
500 : end if
501 : end do
502 :
503 : ! Usually, the q points come in the right order. However, this is not always the case...
504 21 : if(found==0)then
505 :
506 : ! If the EIG2 database file has to be read again, close it, then search for the right q point,
507 : ! from the beginning of the file
508 0 : call ddb_hdr%close()
509 :
510 0 : call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self, msym=msym)
511 :
512 : ! And examine again the EIG2 file. Still, not beyond the previously examined value.
513 0 : found=0
514 0 : do iqpt2=1,iqpt2_previous
515 0 : call ddb_eig2%read_d2eig(ddb_hdr, iqpt2, iqpt2)
516 0 : diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt)
517 0 : if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then
518 : found=1
519 : exit
520 : end if
521 : end do
522 :
523 0 : if(found==0)then
524 : write(message,'(a,3es16.6,2a)')&
525 0 : & 'Was unable to find in EIG2 DDB file, the blok for point ',spqpt(:,iqpt),ch10,&
526 0 : & 'Action: compute the contribution from this point, and merge it in your EIG2 DDB file.'
527 0 : ABI_ERROR(message)
528 : end if
529 :
530 0 : call ddb_hdr%free()
531 :
532 : end if
533 :
534 : ! Put eig2dval in cartesian coordinates
535 : call carttransf(ddb_eig2%flg,ddb_eig2%eig2dval,carflg_eig2,gprimd,iqpt,mband,&
536 21 : & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd)
537 :
538 148 : do imod=1,3*natom
539 :
540 : ! Calculate the derivative
541 7218 : deigr(:,:) = zero
542 7218 : deigi(:,:) = zero
543 7218 : dwtermr(:,:)=zero
544 7218 : dwtermi(:,:)=zero
545 480 : index=0
546 480 : do iatom1=1,natom
547 1560 : do idir1=1,3
548 6984 : do iatom2=1,natom
549 : ! Compute factor for SE term
550 5544 : if(phfreq(imod,iqpt)<tol6)then
551 : factr = zero
552 : else
553 4755 : factr=one/sqrt(amu(typat(iatom1))*amu(typat(iatom2)))/phfreq(imod,iqpt)/amu_emass
554 : end if
555 :
556 23256 : do idir2=1,3
557 16632 : index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1)))
558 :
559 : ! Compute products of polarization vectors
560 : vecr = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)+&
561 16632 : & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod)
562 : veci = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)-&
563 16632 : & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod)
564 :
565 : vec1r = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)+&
566 16632 : & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod)
567 : vec1i = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)-&
568 16632 : & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod)
569 :
570 : vec2r = eigvec(1,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)+&
571 16632 : & eigvec(2,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod)
572 : vec2i = eigvec(2,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)-&
573 16632 : & eigvec(1,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod)
574 :
575 : ! Compute factor for DW term
576 16632 : if(phfreq(imod,iqpt)<tol6)then
577 : fact2r = zero
578 : fact2i = zero
579 : else
580 : fact2r = -wghtq(iqpt)*(vec1r/amu(typat(iatom1)) + vec2r/amu(typat(iatom2)))/phfreq(imod,iqpt)/&
581 14265 : & amu_emass/2 !/norm(idir1)/norm(idir2)
582 : fact2i = -wghtq(iqpt)*(vec1i/amu(typat(iatom1)) + vec2i/amu(typat(iatom2)))/phfreq(imod,iqpt)/&
583 14265 : & amu_emass/2 !/norm(idir1)/norm(idir2)
584 : end if
585 :
586 673920 : multr(:,:) =(ddb_eig2%eig2dval(1,index,:,:)*vecr - ddb_eig2%eig2dval(2,index,:,:)*veci) !/(norm(idir1)*norm(idir2))
587 673920 : multi(:,:) =(ddb_eig2%eig2dval(1,index,:,:)*veci + ddb_eig2%eig2dval(2,index,:,:)*vecr) !/(norm(idir1)*norm(idir2))
588 :
589 :
590 : ! Debye-Waller Term
591 16632 : if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
592 673920 : dwtermr(1:mband,1:nkpt)=dwtermr(1:mband,1:nkpt)+fact2r*eig2dGamma(1,index,:,:)-fact2i*eig2dGamma(2,index,:,:)
593 673920 : dwtermi(1:mband,1:nkpt)=dwtermi(1:mband,1:nkpt)+fact2r*eig2dGamma(2,index,:,:)+fact2i*eig2dGamma(1,index,:,:)
594 : end if
595 :
596 : ! Self-energy Term (Fan)
597 673920 : deigr(1:mband,1:nkpt) = deigr(1:mband,1:nkpt) + wghtq(iqpt)*factr*multr(1:mband,1:nkpt)
598 679464 : deigi(1:mband,1:nkpt) = deigi(1:mband,1:nkpt) + wghtq(iqpt)*factr*multi(1:mband,1:nkpt)
599 :
600 : end do !idir2
601 : end do !iatom2
602 : end do !idir1
603 : end do !iatom1
604 : ! Eigenvalue derivative or broadening
605 141 : if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
606 7218 : dednr(1:mband,1:nkpt,imod,iqpt) = deigr(1:mband,1:nkpt) + dwtermr(1:mband,1:nkpt)
607 7218 : dedni(1:mband,1:nkpt,imod,iqpt) = deigi(1:mband,1:nkpt) + dwtermi(1:mband,1:nkpt)
608 : else if(thmflag==4 .or. thmflag==6 .or. thmflag==8) then
609 0 : dednr(1:mband,1:nkpt,imod,iqpt) = pi*deigr(1:mband,1:nkpt)
610 0 : dedni(1:mband,1:nkpt,imod,iqpt) = pi*deigi(1:mband,1:nkpt)
611 : end if
612 :
613 : end do ! imod
614 : end do !iqpt
615 :
616 7 : call ddb_hdr%close()
617 :
618 :
619 : !=============================================================================
620 : !3) Evaluation of the Eliashberg type spectral function
621 : !and phonon DOS via gaussian broadning
622 : !=============================================================================
623 :
624 7 : if(telphint==1)then
625 7 : ng2f = 500 ! number of frequencies
626 7 : omega_min=zero
627 7 : omega_max=zero
628 28 : do iqpt=1,nqpt
629 148 : do imod=1,3*natom
630 120 : omega_min = min(omega_min,phfreq(imod,iqpt))
631 141 : omega_max = max(omega_max,phfreq(imod,iqpt))
632 : end do
633 : end do
634 :
635 7 : ABI_MALLOC(dos_phon,(ng2f))
636 35 : ABI_MALLOC(g2f,(mband,nkpt,ng2f))
637 28 : ABI_MALLOC(tmpg2f,(mband,nkpt,ng2f))
638 7 : ABI_MALLOC(tmpphondos,(ng2f))
639 :
640 7 : write(std_out,'(a,es13.6)') 'omega_min :', omega_min
641 7 : write(std_out,'(a,es13.6)') 'omega_max :', omega_max
642 7 : write(std_out,'(a,i8)') 'ng2f :', ng2f
643 :
644 7 : omega_max = omega_max + 0.1 * omega_max
645 7 : domega = (omega_max-omega_min)/(ng2f-one)
646 :
647 7 : gaussprefactor = sqrt(piinv) / g2fsmear
648 7 : gaussfactor = one / g2fsmear
649 :
650 227507 : g2f(:,:,:) = zero
651 3507 : dos_phon(:) = zero
652 :
653 28 : do iqpt=1,nqpt
654 148 : do imod=1,3*natom
655 60120 : omega = omega_min
656 3609120 : tmpg2f(:,:,:) = zero
657 60120 : tmpphondos(:) = zero
658 60120 : do iomega=1,ng2f
659 60000 : xx = (omega-phfreq(imod,iqpt))*gaussfactor
660 60000 : gaussval = gaussprefactor*exp(-xx*xx)
661 3609000 : tmpg2f(:,:,iomega) = tmpg2f(:,:,iomega) + gaussval*dednr(:,:,imod,iqpt)
662 60000 : tmpphondos(iomega) = tmpphondos(iomega) + gaussval
663 60120 : omega = omega+domega
664 : end do
665 :
666 3609120 : g2f(:,:,:) = g2f(:,:,:) + tmpg2f(:,:,:)
667 60141 : dos_phon(:) = dos_phon(:) + tmpphondos(:)
668 :
669 : end do !imod
670 : end do !iqpt
671 :
672 3507 : dos_phon(:) = dos_phon(:) / nqpt
673 :
674 : ! output the g2f
675 255 : kpnt(:,:,1) = ddb_eig2%kpt(:,:)
676 7 : unit_g2f = 108
677 7 : call outg2f(domega,omega_min,omega_max,elph_base_name,g2f,g2fsmear,kpnt,mband,ng2f,nkpt,nqpt,1,telphint,unit_g2f)
678 :
679 : ! output the phonon DOS
680 7 : unit_phdos = 108
681 7 : call outphdos(domega,dos_phon,omega_min,omega_max,elph_base_name,g2fsmear,ng2f,nqpt,1,telphint,unit_g2f)
682 :
683 :
684 7 : ABI_FREE(dos_phon)
685 7 : ABI_FREE(g2f)
686 7 : ABI_FREE(tmpg2f)
687 7 : ABI_FREE(tmpphondos)
688 :
689 : end if !telphint
690 :
691 : !=======================================================================
692 : !4) Evaluation of the Eliashberg type spectral function
693 : !and phonon DOS via improved tetrahedron method
694 : !=======================================================================
695 :
696 7 : if(telphint==0)then
697 :
698 : ! make dimension-ful rprimd and gprimd for transformation of derivatives to cartesian coordinates.
699 0 : call mkrdim(acell,rprim,rprimd)
700 0 : call matr3inv(rprimd,gprimd)
701 :
702 : ! Q point Grid
703 0 : qpt_full(:,:) = ddb%qpt(1:3,:)
704 :
705 : ! Trivial Q point index
706 0 : do iqpt=1,nqpt
707 0 : indqpt(iqpt)=iqpt
708 0 : qptnrm(iqpt)= qpt_full(1,iqpt)*qpt_full(1,iqpt)+qpt_full(2,iqpt)*qpt_full(2,iqpt)+qpt_full(3,iqpt)*qpt_full(3,iqpt)
709 : end do
710 :
711 : ! Build qlatt from scratch (for 5.7)
712 0 : tol = 0.1_dp
713 0 : ilatt = 0
714 0 : call sort_dp(nqpt,qptnrm,indqpt,tol)
715 :
716 0 : do iqpt1=1,nqpt-2
717 0 : mesh(1:3,1) = qpt_full(1:3,indqpt(iqpt1))
718 0 : do iqpt2=iqpt1+1,nqpt-1
719 0 : mesh(1:3,2)= qpt_full(1:3,indqpt(iqpt2))
720 0 : do iqpt3=iqpt2+1,nqpt
721 0 : mesh(1:3,3)= qpt_full(1:3,indqpt(iqpt3))
722 : det = mesh(1,1)*mesh(2,2)*mesh(3,3) + mesh(1,2)*mesh(2,3)*mesh(3,1) + mesh(1,3)*mesh(2,1)*mesh(3,2) &
723 0 : & -mesh(3,1)*mesh(2,2)*mesh(1,3) - mesh(3,2)*mesh(2,3)*mesh(1,1) - mesh(3,3)*mesh(2,1)*mesh(1,2)
724 0 : invdet = one/det
725 0 : if (abs(nint(invdet))==nqpt .and. abs(invdet)-nqpt < tol) then
726 0 : ilatt = 1
727 0 : qlatt(:,:) = mesh(:,:)
728 : exit
729 : end if
730 : end do
731 0 : if(ilatt==1) exit
732 : end do
733 0 : if(ilatt==1) exit
734 : end do
735 :
736 : ! error message if qlatt not found and stop
737 0 : if(ilatt==0) then
738 : write(message, '(a,a)' ) &
739 0 : & ' Could not find homogeneous basis vectors for Q point grid ',ch10
740 0 : call wrtout(std_out,message,'COLL')
741 0 : call wrtout(ab_out,message,'COLL')
742 0 : ABI_ERROR("Aborting now")
743 : end if
744 :
745 : ! test if qlatt is righthanded and possibly fixe it
746 0 : if(invdet < 0) then
747 0 : tempqlatt(:) = qlatt(:,2)
748 0 : qlatt(:,2) = qlatt(:,1)
749 0 : qlatt(:,1) = tempqlatt(:)
750 : end if
751 :
752 0 : write(std_out,*) 'qlatt',qlatt
753 :
754 : ! test if qlatt generates all Q points TO DO
755 :
756 : ! Get tetrahedra, ie indexes of the full kpoints at their summits
757 0 : call init_tetra(indqpt,gprimd,qlatt,qpt_full,nqpt, tetrahedra, ierr, errstr, xmpi_comm_self)
758 : !call htetra_init(tetra, indqpt, gprimd, qlatt, qpt_full, nqpt, kpt_ibz, nkpt_ibz, ierr, errstr, xmpi_comm_self
759 0 : ABI_CHECK(ierr==0,errstr)
760 :
761 : rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
762 : & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
763 0 : & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
764 :
765 : ! Calculate weights for phonon DOS
766 : ! Special precautions must be taking for Gamma point
767 : ! because of non-analytic term.
768 : ! Non-analyticity must be taken out and treated separatly.
769 :
770 0 : nene = 100 !nene=number of energies for DOS
771 0 : enemin = minval(phfreq)
772 0 : enemax = maxval(phfreq)
773 0 : deltaene = (enemax-enemin)/dble(nene-1)
774 : ! redefine enemin enemax to be at rounded multiples of deltaene
775 : ! enemin = elph_ds%fermie - dble(ifermi)*deltaene
776 : ! enemax = elph_ds%fermie + dble(nene-ifermi-1)*deltaene
777 :
778 0 : ABI_MALLOC(tweight,(nqpt,nene))
779 0 : ABI_MALLOC(dtweightde,(nqpt,nene))
780 0 : ABI_MALLOC(intweight,(3*natom,nqpt,nene))
781 0 : ABI_MALLOC(indtweightde,(3*natom,nqpt,nene))
782 :
783 0 : do iband=1,3*natom
784 0 : eigen_in(:) = phfreq(iband,:)
785 :
786 : ! calculate general integration weights at each irred kpoint
787 : ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
788 : call get_tetra_weight(eigen_in,enemin,enemax,&
789 : & one,nene,nqpt,tetrahedra,bcorr0,&
790 0 : & tweight,dtweightde,xmpi_comm_self)
791 :
792 0 : intweight(iband,:,:) = tweight(:,:)
793 0 : indtweightde(iband,:,:) = dtweightde(:,:)
794 :
795 : end do !iband
796 :
797 : ! intdtweightse(nband,nqpt,nene) represents the weight in each energy bin for every kpt and every band
798 : ! So phonon DOS is calculated (neglecting the non-analyticity contribution for now !!!)
799 :
800 0 : ABI_MALLOC(total_dos,(nene))
801 0 : ABI_MALLOC(g2f,(mband,nkpt,nene))
802 :
803 0 : total_dos(:) = zero
804 0 : do iband=1,3*natom
805 0 : do iqpt=1,nqpt
806 0 : total_dos(:) = total_dos + indtweightde(iband,iqpt,:)
807 : end do
808 : end do
809 :
810 : ! For the g2f function
811 : ! Right now for one electronic band and one K point: dednr(1:mband,1:nkpt,imod,iqpt)
812 : ! Once again must pay close attention to the Gamma point
813 0 : g2f(:,:,:) = zero
814 0 : do ii=1,mband
815 0 : do ikpt=1,nkpt
816 0 : do iband=1,3*natom
817 0 : do iqpt=1,nqpt
818 0 : g2f(ii,ikpt,:) = g2f(ii,ikpt,:) + dednr(ii,ikpt,iband,iqpt) * indtweightde(iband,iqpt,:)
819 : end do
820 : end do
821 : end do
822 : end do
823 :
824 : ! output the g2f
825 0 : unit_g2f = 108
826 0 : call outg2f(deltaene,enemin,enemax,elph_base_name,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,tetrahedra%ntetra,telphint,unit_g2f)
827 :
828 : ! output the phonon DOS
829 0 : unit_phdos = 108
830 0 : call outphdos(deltaene,total_dos,enemin,enemax,elph_base_name,g2fsmear,nene,nqpt,tetrahedra%ntetra,telphint,unit_g2f)
831 :
832 0 : ABI_FREE(tweight)
833 0 : ABI_FREE(dtweightde)
834 0 : ABI_FREE(intweight)
835 0 : ABI_FREE(indtweightde)
836 0 : ABI_FREE(total_dos)
837 0 : ABI_FREE(g2f)
838 : end if !telphint
839 :
840 : !=======================================================================
841 : !5) direct evaluation of thermal corrections
842 : !=======================================================================
843 :
844 : !open TBS file
845 7 : outfile = trim(elph_base_name)//"_TBS"
846 7 : if (open_file(outfile,message,newunit=unitout,form='formatted',status='unknown') /= 0) then
847 0 : ABI_ERROR(message)
848 : end if
849 7 : write(unitout,'(a)')'thmeig: Thermal Eigenvalue corrections (eV)'
850 :
851 1227 : slope(:,:,:) = zero
852 1227 : zeropoint(:,:,:) = zero
853 : !Loop on temperatures
854 77 : do itemper= 1, ntemper
855 70 : tmp=tempermin+temperinc*float(itemper-1)
856 12270 : thmeigen(:,:,:) = zero
857 :
858 : ! Sum on all phonon wavevectors and modes
859 280 : do iqpt=1,nqpt
860 1480 : do imod=1,3*natom
861 :
862 : ! Bose-Einstein distribution
863 : ! jmb overflow with exp(). So, select bosein to be still significant wrt half
864 1200 : if(phfreq(imod,iqpt)<tol6 .or. (phfreq(imod,iqpt)/(kb_HaK*tmp)) > -log(tol16))then
865 : bosein = zero
866 : else
867 969 : bosein = one/(exp(phfreq(imod,iqpt)/(kb_HaK*tmp))-one)
868 : end if
869 :
870 : ! Calculate total
871 72180 : thmeigen(1,1:mband,1:nkpt) = thmeigen(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(bosein+half)
872 72180 : thmeigen(2,1:mband,1:nkpt) = thmeigen(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(bosein+half)
873 :
874 1410 : if(itemper==1)then
875 : ! Calculate slope of linear regime
876 120 : if(phfreq(imod,iqpt)<tol6)then
877 : slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt)
878 : slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt)
879 : else
880 5845 : slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt))
881 5845 : slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt))
882 : end if
883 : ! Calculate zero-point renormalization
884 7218 : zeropoint(1,1:mband,1:nkpt) = zeropoint(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*half
885 7218 : zeropoint(2,1:mband,1:nkpt) = zeropoint(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*half
886 :
887 : end if
888 : end do ! imod
889 : end do !iqpt
890 :
891 : ! Write temperature independent results
892 70 : if(itemper==1)then
893 7 : write(unitout,'(a)')'Temperature independent results (zero-point renormalization and slope)'
894 69 : do ikpt=1,nkpt
895 62 : write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
896 455 : do iband=1,mband
897 386 : write(unitout,'(4d22.14)') Ha_eV*zeropoint(1,iband,ikpt),Ha_eV*zeropoint(2,iband,ikpt),&
898 834 : & Ha_eV*slope(1,iband,ikpt),Ha_eV*slope(2,iband,ikpt)
899 : end do
900 : end do
901 7 : write(unitout,'(a)')'Temperature dependent corrections'
902 : end if
903 : ! Write result in file for each temperature
904 70 : write(unitout,'(a,es10.3,a)')'T :', tmp,' K'
905 697 : do ikpt=1,nkpt
906 620 : write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
907 4550 : do iband=1,mband
908 4480 : write(unitout,'(2d22.14)') Ha_eV*thmeigen(1,iband,ikpt), Ha_eV*thmeigen(2,iband,ikpt)
909 : end do
910 : end do
911 : end do !itemper
912 :
913 7 : close(unitout)
914 :
915 : !Write temperature-independent results to the main output file
916 7 : write(iout,'(a)')' '
917 567 : write(iout,'(80a)') ('-',ii=1,80)
918 7 : write(iout,'(a)')' '
919 7 : write(iout,'(a)')' Electron-phonon change of electronic structure.'
920 7 : write(iout,'(a)')' The temperature-dependent values are written in the _TBS file.'
921 7 : write(iout,'(a)')' Here follows, for each electronic wavevector and band :'
922 7 : write(iout,'(a)')' zero-point renormalisation (Ha) and linear slope (Ha/Kelvin)'
923 69 : do ikpt=1,nkpt
924 62 : write(iout,'(2a,i6,a,3es16.6)')ch10,' Kpt number ',ikpt,', with reduced coordinates :',kpnt(:,ikpt,1)
925 455 : do iband=1,mband
926 448 : write(iout,'(i6,2es20.6)') iband,zeropoint(1,iband,ikpt),slope(1,iband,ikpt)
927 : end do
928 : end do
929 :
930 7 : ABI_FREE(typat)
931 7 : ABI_FREE(zion)
932 7 : ABI_FREE(amu)
933 7 : ABI_FREE(xcart)
934 7 : ABI_FREE(xred)
935 7 : ABI_FREE(symafm)
936 7 : ABI_FREE(spinat)
937 7 : ABI_FREE(symrel)
938 7 : ABI_FREE(symrec)
939 7 : ABI_FREE(indsym)
940 7 : ABI_FREE(tnons)
941 7 : ABI_FREE(deigi)
942 7 : ABI_FREE(deigr)
943 7 : ABI_FREE(dwtermi)
944 7 : ABI_FREE(dwtermr)
945 7 : ABI_FREE(multi)
946 7 : ABI_FREE(multr)
947 7 : ABI_FREE(slope)
948 7 : ABI_FREE(thmeigen)
949 7 : ABI_FREE(zeropoint)
950 :
951 7 : ABI_FREE(dedni)
952 7 : ABI_FREE(dednr)
953 7 : if(allocated(indqpt)) then
954 7 : ABI_FREE(indqpt)
955 : end if
956 7 : ABI_FREE(eigen_in)
957 7 : ABI_FREE(qpt_full)
958 7 : ABI_FREE(qptnrm)
959 7 : ABI_FREE(wghtq)
960 7 : ABI_FREE(spqpt)
961 7 : ABI_FREE(eigvec)
962 7 : ABI_FREE(phfreq)
963 :
964 7 : ABI_FREE(eig2dGamma)
965 7 : ABI_FREE(kpnt)
966 7 : ABI_FREE(carflg_eig2)
967 :
968 7 : call ddb_eig2%free()
969 7 : call destroy_tetra(tetrahedra)
970 :
971 7 : end subroutine thmeig
972 : !!***
973 :
974 : !!****f* m_thmeig/outphdos
975 : !! NAME
976 : !! outphdos
977 : !!
978 : !! FUNCTION
979 : !! Print out phonon density of states
980 : !!
981 : !! INPUTS
982 : !! deltaene = step on energy/frequency grid, in Hartree
983 : !! dos_phon = phonon DOS calculated on a grid
984 : !! enemin = minimal frequency
985 : !! enemax = maximal frequency
986 : !! filnam = file name for output to disk
987 : !! g2fsmear = smearing width
988 : !! nene = number of points on energy axis
989 : !! nqpt = number of q-points
990 : !! ntetra = number of tetrahedra, if tetrahedron interpolation is used
991 : !! telphint = flag for el-phonon interpolation method (to indicate Gaussian or tetrahedron integration)
992 : !! unit_phdos = unit for phonon DOS output
993 : !!
994 : !!
995 : !! OUTPUT
996 : !! only write
997 : !!
998 : !! SIDE EFFECTS
999 : !!
1000 : !! NOTES
1001 : !! FIXME
1002 : !! overcomplete inputs. Eliminate unit_phdos (just filnam) and deltaene (gotten from max-min/nene)
1003 : !!
1004 : !! SOURCE
1005 :
1006 7 : subroutine outphdos(deltaene,dos_phon,enemin,enemax,filnam,g2fsmear,nene,nqpt,ntetra,telphint,unit_phdos)
1007 :
1008 : !Arguments ------------------------------------
1009 : !scalars
1010 : integer,intent(in) :: nene,nqpt,ntetra,telphint,unit_phdos
1011 : character(len=fnlen),intent(in) :: filnam
1012 : real(dp) :: deltaene,enemin,enemax,g2fsmear
1013 : !arrays
1014 : real(dp) :: dos_phon(nene)
1015 :
1016 : !Local variables-------------------------------
1017 : !scalars
1018 : integer :: iomega,iost,step10
1019 : real(dp) :: dos_effective,omega
1020 : character(len=fnlen) :: outfile
1021 : character(len=500) :: message
1022 : !arrays
1023 :
1024 : ! *************************************************************************
1025 :
1026 7 : outfile = trim(filnam) // '_PDS'
1027 7 : write(message, '(3a)')ch10,&
1028 14 : & ' Will write phonon DOS in file ',trim(outfile)
1029 7 : call wrtout(ab_out,message,'COLL')
1030 7 : call wrtout(std_out,message,'COLL')
1031 :
1032 7 : write(message, '(4a)')ch10,&
1033 7 : & ' For checking purposes, write ten values in the present file.',ch10,&
1034 14 : & ' Index Energy (in Ha) DOS '
1035 7 : call wrtout(ab_out,message,'COLL')
1036 7 : call wrtout(std_out,message,'COLL')
1037 :
1038 7 : open (unit=unit_phdos,file=outfile,status='replace',iostat=iost)
1039 7 : if (iost /= 0) then
1040 0 : write (message,'(3a)')' Opening file ',trim(outfile),' as new'
1041 0 : ABI_ERROR(message)
1042 : end if
1043 :
1044 7 : write (unit_phdos,'(a)') '#'
1045 7 : write (unit_phdos,'(a)') '# ABINIT package : phonon DOS file'
1046 7 : write (unit_phdos,'(a)') '#'
1047 7 : write (unit_phdos,'(a,i10)') '# Number of Qpoints integrated over : ', nqpt
1048 7 : write (unit_phdos,'(a,i10)') '# Number of energy points : ', nene
1049 7 : write (unit_phdos,'(a,es16.6,a,es16.6,a)') '# between omega_min = ', enemin, &
1050 14 : & ' Ha and omega_max = ', enemax, ' Ha'
1051 7 : if(telphint==1)then
1052 7 : write (unit_phdos,'(a,es16.6)') '# The smearing width for gaussians is ', g2fsmear
1053 : end if
1054 7 : if(telphint==0)then
1055 0 : write (unit_phdos,'(a,i10)') '# Number of tetrahedrons', ntetra
1056 : end if
1057 7 : write (unit_phdos,'(a)') '#'
1058 7 : write (unit_phdos,'(a)') '# Index Energy (in Ha) DOS '
1059 :
1060 7 : omega = enemin
1061 3507 : do iomega=1,nene
1062 3500 : dos_effective=dos_phon(iomega)
1063 3500 : if(abs(dos_effective)<tol16)then
1064 2738 : dos_effective=zero
1065 : end if
1066 3500 : step10=nene/10
1067 3500 : if(mod(iomega,step10)==1)write (std_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective
1068 3500 : if(mod(iomega,step10)==1)write (ab_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective
1069 3500 : write (unit_phdos, '(i10,es18.6,es18.6)')iomega, omega, dos_effective
1070 3507 : omega=omega+deltaene
1071 : end do
1072 :
1073 7 : close (unit=unit_phdos)
1074 :
1075 7 : end subroutine outphdos
1076 : !!***
1077 :
1078 : !!****f* m_thmeig/outg2f
1079 : !! NAME
1080 : !! outg2f
1081 : !!
1082 : !! FUNCTION
1083 : !! Output g2f function to file. FIXME: Paul, please explain what g2f is.
1084 : !! Probably a variant on the Eliashberg spectral function a2F
1085 : !!
1086 : !! INPUTS
1087 : !!
1088 : !! OUTPUT
1089 : !! only write
1090 : !!
1091 : !! SOURCE
1092 :
1093 7 : subroutine outg2f(deltaene,enemin,enemax,filnam,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f)
1094 :
1095 : !Arguments ------------------------------------
1096 : !scalars
1097 : integer,intent(in) :: mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f
1098 : character(len=fnlen),intent(in) :: filnam
1099 : real(dp) :: deltaene,enemin,enemax,g2fsmear
1100 : !arrays
1101 : real(dp) :: g2f(mband,nkpt,nene),kpnt(3,nkpt,nqpt)
1102 :
1103 : !Local variables-------------------------------
1104 : !scalars
1105 : integer :: iband,ikpt,iomega,iost
1106 : real(dp) :: omega
1107 : character(len=fnlen) :: outfile
1108 : character(len=500) :: message
1109 : !arrays
1110 :
1111 : ! *************************************************************************
1112 :
1113 : !output the g2f
1114 7 : outfile = trim(filnam) // '_G2F'
1115 7 : open (unit=unit_g2f,file=outfile,status='unknown',iostat=iost)
1116 7 : if (iost /= 0) then
1117 0 : write (message,'(3a)')' thmeig : ERROR- opening file ',trim(outfile),' as new'
1118 0 : ABI_ERROR(message)
1119 : end if
1120 :
1121 7 : write(std_out,*) ' g2f function'
1122 7 : write (unit_g2f,'(a)') '#'
1123 7 : write (unit_g2f,'(a)') '# ABINIT package : g2f file'
1124 7 : write (unit_g2f,'(a)') '#'
1125 7 : write (unit_g2f,'(a,I10)') '# number of qpoints integrated over : ', nqpt
1126 7 : write (unit_g2f,'(a,I10)') '# number of energy points : ', nene
1127 7 : write (unit_g2f,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', enemin, &
1128 14 : & ' Ha and omega_max = ', enemax, ' Ha'
1129 7 : if(telphint==1)then
1130 7 : write (unit_g2f,'(a,E16.6)') '# and the smearing width for gaussians is ', g2fsmear
1131 7 : write (unit_g2f,'(a)') '#'
1132 : end if
1133 7 : if(telphint==0)then
1134 0 : write (unit_g2f,'(a,I10)') '# number of tetrahedrons', ntetra
1135 0 : write (unit_g2f,'(a)') '#'
1136 : end if
1137 :
1138 : !Write only the a2f function for the first K point
1139 : !ikpt=1
1140 69 : do ikpt=1,nkpt
1141 62 : write(unit_g2f,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
1142 455 : do iband=1,mband
1143 386 : write(unit_g2f,*) 'band :', iband
1144 386 : omega = enemin
1145 193448 : do iomega=1,nene
1146 193000 : write (unit_g2f,*) omega*Ha_eV*1000, g2f(iband, ikpt,iomega)
1147 193386 : omega=omega+deltaene
1148 : end do
1149 : end do
1150 : end do
1151 :
1152 7 : close (unit=unit_g2f)
1153 :
1154 7 : end subroutine outg2f
1155 : !!***
1156 :
1157 : end module m_thmeig
1158 : !!***
|
