Branch data 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 : : !!***
|
