TABLE OF CONTENTS


ABINIT/m_frohlichmodel [ Modules ]

[ Top ] [ Modules ]

NAME

  m_frohlichmodel

FUNCTION

  Compute ZPR, temperature-dependent electronic structure, and other properties
  using the Frohlich model 

COPYRIGHT

  Copyright (C) 2018 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_frohlichmodel
29 
30  use defs_basis
31  use m_errors
32  use m_abicore
33 
34  implicit none
35 
36  private
37 
38  public :: frohlichmodel
39 
40 contains

m_frohlichmodel/frohlichmodel [ Functions ]

[ Top ] [ m_frohlichmodel ] [ Functions ]

NAME

  frohlichmodel

FUNCTION

 Main routine to compute properties based on the Frohlich model

INPUTS

 cryst<crystal_t>=Structure defining the unit cell (geometry, atomic positions and symmetry operations in real and reciprocal space)
 dtfil<datafiles_type>=Variables related to files.
 dtset<dataset_type>=All input variables for this dataset.
 ebands<ebands_t>=Electronic band structure information
 efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point
 efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure
   efmasval(:,:)%eig2_diag band curvature double tensor
 ifc<ifc_type>=contains the dynamical matrix and the IFCs.

PARENTS

 eph

CHILDREN

SOURCE

 67 subroutine frohlichmodel(cryst,dtfil,dtset,ebands,efmasdeg,efmasval,ifc)
 68 
 69  use defs_basis
 70  use defs_datatypes
 71  use defs_abitypes
 72  use m_abicore
 73 !use m_xmpi
 74 !use m_xomp
 75  use m_errors
 76 !use m_hdr
 77  use m_crystal
 78  use m_crystal_io
 79  use m_ebands
 80  use m_efmas_defs
 81  use m_ifc
 82 
 83  use m_gaussian_quadrature, only : cgqf
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'frohlichmodel'
 89 !End of the abilint section
 90 
 91  implicit none
 92 
 93 !Arguments ------------------------------------
 94 !scalars
 95  type(crystal_t),intent(in) :: cryst
 96  type(datafiles_type),intent(in) :: dtfil
 97  type(dataset_type),intent(in) :: dtset
 98  type(ebands_t),intent(in) :: ebands
 99  type(ifc_type),intent(in) :: ifc
100 !arrays
101  type(efmasdeg_type), intent(in) :: efmasdeg(:)
102  type(efmasval_type), intent(in) :: efmasval(:,:)
103 
104 !Local variables ------------------------------
105 !scalars
106  logical :: sign_warn
107  integer :: deg_dim,iband,ideg,ikpt,imode,info,iphi,iqdir,itheta
108  integer :: jband,lwork,nphi,nqdir,ntheta
109  real(dp) :: angle_phi,cosph,costh,sinph,sinth,weight,weight_phi
110  real(dp) :: zpr_frohlich,zpr_q0_avg,zpr_q0_fact
111  character(len=500) :: msg
112 !arrays
113  logical, allocatable :: saddle_warn(:), start_eigf3d_pos(:)
114  real(dp) :: kpt(3),unit_r(3)
115  real(dp), allocatable :: eigenval(:), rwork(:), unit_qdir(:,:)
116  real(dp), allocatable :: m_avg(:), m_avg_frohlich(:)
117  real(dp), allocatable :: gq_points_th(:),gq_weights_th(:)
118  real(dp), allocatable :: gq_points_cosph(:),gq_points_sinph(:)
119  real(dp), allocatable :: weight_qdir(:)
120  real(dp), allocatable :: polarity_qdir(:,:,:)
121  real(dp), allocatable :: proj_polarity_qdir(:,:)
122  real(dp), allocatable :: zpr_q0_phononfactor_qdir(:)
123  real(dp), allocatable :: frohlich_phononfactor_qdir(:)
124  real(dp), allocatable :: phfrq_qdir(:,:)
125  real(dp), allocatable :: dielt_qdir(:)
126  real(dp), allocatable :: zpr_frohlich_avg(:)
127  complex(dpc), allocatable :: eigenvec(:,:), work(:)
128  complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:)
129  complex(dpc), allocatable :: f3d(:,:)
130 
131 !************************************************************************
132 
133  !!! Initialization of integrals 
134  ntheta   = dtset%efmas_ntheta
135  nphi     = 2*ntheta
136  nqdir     = nphi*ntheta
137 
138  ABI_ALLOCATE(gq_points_th,(ntheta))
139  ABI_ALLOCATE(gq_weights_th,(ntheta))
140  ABI_ALLOCATE(gq_points_cosph,(nphi))
141  ABI_ALLOCATE(gq_points_sinph,(nphi))
142 
143  ABI_ALLOCATE(unit_qdir,(3,nqdir))
144  ABI_ALLOCATE(weight_qdir,(nqdir))
145 
146  call cgqf(ntheta,1,zero,zero,zero,pi,gq_points_th,gq_weights_th)
147  weight_phi=two*pi/real(nphi,dp)
148  do iphi=1,nphi
149    angle_phi=weight_phi*(iphi-1)
150    gq_points_cosph(iphi)=cos(angle_phi)
151    gq_points_sinph(iphi)=sin(angle_phi)
152  enddo
153  nqdir=0
154  do itheta=1,ntheta
155    costh=cos(gq_points_th(itheta))
156    sinth=sin(gq_points_th(itheta))
157    weight=gq_weights_th(itheta)*weight_phi*sinth
158    do iphi=1,nphi
159      cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
160      nqdir=nqdir+1
161 
162      unit_qdir(1,nqdir)=sinth*cosph
163      unit_qdir(2,nqdir)=sinth*sinph
164      unit_qdir(3,nqdir)=costh
165      weight_qdir(nqdir)=weight
166 
167    enddo
168  enddo
169 
170  ABI_DEALLOCATE(gq_points_th)
171  ABI_DEALLOCATE(gq_weights_th)
172  ABI_DEALLOCATE(gq_points_cosph)
173  ABI_DEALLOCATE(gq_points_sinph)
174 
175  ABI_ALLOCATE(polarity_qdir,(3,3*cryst%natom,nqdir))
176  ABI_ALLOCATE(proj_polarity_qdir,(3*cryst%natom,nqdir))
177  ABI_ALLOCATE(zpr_q0_phononfactor_qdir,(nqdir))
178  ABI_ALLOCATE(frohlich_phononfactor_qdir,(nqdir))
179  ABI_ALLOCATE(phfrq_qdir,(3*cryst%natom,nqdir))
180  ABI_ALLOCATE(dielt_qdir,(nqdir))
181 
182  !Compute phonon frequencies and mode-polarity for each qdir
183  call ifc_calcnwrite_nana_terms(ifc, cryst, nqdir, unit_qdir, &
184 &  phfrq2l=phfrq_qdir, polarity2l=polarity_qdir)
185 
186  !Compute dielectric tensor for each qdir
187  do iqdir=1,nqdir
188    dielt_qdir(iqdir)=DOT_PRODUCT(unit_qdir(:,iqdir),MATMUL(ifc%dielt(:,:),unit_qdir(:,iqdir)))
189  enddo
190 
191  !Compute projection of mode-polarity on qdir, and other derived quantities summed over phonon branches for each iqdir.
192  !Note that acoustic modes are discarded (imode sum starts only from 4)
193  zpr_q0_phononfactor_qdir=zero
194  zpr_q0_avg=zero
195  frohlich_phononfactor_qdir=zero
196  do iqdir=1,nqdir
197    do imode=4,3*cryst%natom
198      proj_polarity_qdir(imode,iqdir)=DOT_PRODUCT(unit_qdir(:,iqdir),polarity_qdir(:,imode,iqdir))
199      zpr_q0_phononfactor_qdir(iqdir)=zpr_q0_phononfactor_qdir(iqdir)+&
200 &      proj_polarity_qdir(imode,iqdir)**2 / phfrq_qdir(imode,iqdir) **2
201      frohlich_phononfactor_qdir(iqdir)=frohlich_phononfactor_qdir(iqdir)+&
202 &      proj_polarity_qdir(imode,iqdir)**2 / phfrq_qdir(imode,iqdir) **(three*half)
203    enddo
204    zpr_q0_avg=zpr_q0_avg+&
205 &    weight_qdir(iqdir)*zpr_q0_phononfactor_qdir(iqdir)/dielt_qdir(iqdir)**2
206  enddo
207  zpr_q0_avg=zpr_q0_avg*quarter*piinv
208  zpr_q0_fact=zpr_q0_avg*two*pi*pi*(three*quarter*piinv)**third*cryst%ucvol**(-four*third)
209 
210 !DEBUG
211 ! do iqdir=1,nqdir,513
212 !   write(std_out,'(a,3f8.4,3es12.4)')' unit_qdir,dielt_qdir,zpr_q0_phononfactor_qdir,frohlich_phononfactor=',&
213 !&    unit_qdir(:,iqdir),dielt_qdir(iqdir),zpr_q0_phononfactor_qdir(iqdir),frohlich_phononfactor_qdir(iqdir)
214 !   do imode=1,3*cryst%natom
215 !     write(std_out,'(a,i5,6es12.4)')'   imode,phfrq_qdir,phfrq(cmm1),polarity_qdir=',&
216 !&     imode,phfrq_qdir(imode,iqdir),phfrq_qdir(imode,iqdir)*Ha_cmm1,polarity_qdir(:,imode,iqdir),proj_polarity_qdir(imode,iqdir)
217 !   enddo
218 ! enddo
219 ! write(std_out,'(2a,3es12.4)')ch10,&
220 !& ' zpr_q0_avg, zpr_q0_fact, zpr_q0_fact (eV) =',zpr_q0_avg, zpr_q0_fact, zpr_q0_fact*Ha_eV
221 !ENDDEBUG
222 
223  write(ab_out,'(6a,f14.6,a,f14.6,a)') ch10,&
224 &  ' Rough correction to the ZPR, to take into account the missing q=0 piece using Frohlich model:',ch10,&
225 &  ' (+ for occupied states, - for unoccupied states) * zpr_q0_fact / (Nqpt_full_bz)**(1/3) ',ch10,&
226 &  ' where Nqpt_full_bz=number of q wavevectors in full BZ, and zpr_q0_fact=',zpr_q0_fact,' Ha=',zpr_q0_fact*Ha_eV,' eV'
227 
228  !Compute effective masses, and integrate the Frohlich model
229  do ikpt=1,dtset%nkpt
230 
231    kpt(:)=dtset%kptns(:,ikpt)
232    do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
233 
234      deg_dim    = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
235 
236      ABI_ALLOCATE(eig2_diag_cart,(3,3,deg_dim,deg_dim))
237 
238      !Convert eig2_diag to cartesian coordinates
239      do iband=1,deg_dim
240         do jband=1,deg_dim
241           eig2_diag_cart(:,:,iband,jband)=efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)
242           eig2_diag_cart(:,:,iband,jband)=&
243 &           matmul(matmul(cryst%rprimd,eig2_diag_cart(:,:,iband,jband)),transpose(cryst%rprimd))/two_pi**2
244         enddo
245      enddo
246 
247      ABI_ALLOCATE(f3d,(deg_dim,deg_dim))
248      ABI_ALLOCATE(m_avg,(deg_dim))
249      ABI_ALLOCATE(m_avg_frohlich,(deg_dim))
250      ABI_ALLOCATE(zpr_frohlich_avg,(deg_dim))
251      ABI_ALLOCATE(eigenval,(deg_dim))
252      ABI_ALLOCATE(saddle_warn,(deg_dim))
253      ABI_ALLOCATE(start_eigf3d_pos,(deg_dim))
254 
255      m_avg=zero
256      m_avg_frohlich=zero
257      saddle_warn=.false.
258 
259      !Initializations for the diagonalization routine 
260      if(deg_dim>1)then
261 
262        ABI_ALLOCATE(eigenvec,(deg_dim,deg_dim))
263        lwork=-1
264        ABI_ALLOCATE(rwork,(3*deg_dim-2))
265        ABI_ALLOCATE(work,(1))
266        call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
267        lwork=int(work(1))
268        ABI_DEALLOCATE(work)
269        ABI_ALLOCATE(work,(lwork))
270 
271      endif
272 
273      !Perform the integral over the sphere
274      zpr_frohlich_avg=zero
275      do iqdir=1,nqdir
276        do iband=1,deg_dim
277          do jband=1,deg_dim
278            f3d(iband,jband)=DOT_PRODUCT(unit_qdir(:,iqdir),MATMUL(eig2_diag_cart(:,:,iband,jband),unit_qdir(:,iqdir)))
279          enddo
280        enddo
281  
282        if(deg_dim==1)then
283          eigenval(1)=f3d(1,1)
284        else
285          eigenvec = f3d ; eigenval = zero
286          work=zero      ; rwork=zero
287          call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
288        endif
289 
290        m_avg = m_avg + weight_qdir(iqdir)*eigenval
291        m_avg_frohlich = m_avg_frohlich + weight_qdir(iqdir)/(abs(eigenval))**half
292        zpr_frohlich_avg = zpr_frohlich_avg + &
293 &        weight_qdir(iqdir) * frohlich_phononfactor_qdir(iqdir)/((abs(eigenval))**half *dielt_qdir(iqdir)**2)
294 
295        if(iqdir==1) start_eigf3d_pos = (eigenval > 0)
296 
297        do iband=1,deg_dim
298          if(start_eigf3d_pos(iband) .neqv. (eigenval(iband)>0)) then
299            saddle_warn(iband)=.true.
300          end if
301        end do
302 
303      enddo
304 
305      if(deg_dim>1)then
306        ABI_DEALLOCATE(eigenvec)
307        ABI_DEALLOCATE(rwork)
308        ABI_DEALLOCATE(work)
309      endif
310 
311      m_avg = quarter*piinv*m_avg
312      m_avg = one/m_avg
313 
314      m_avg_frohlich = quarter*piinv * m_avg_frohlich
315      m_avg_frohlich = m_avg_frohlich**2
316 
317      zpr_frohlich_avg = quarter*piinv * zpr_frohlich_avg
318 
319      if(deg_dim==1)then
320        write(ab_out,'(2a,3(f6.3,a),i5)')ch10,&
321 &        ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), band ',&
322 &        efmasdeg(ikpt)%degs_bounds(1,ideg)
323      else
324        write(ab_out,'(2a,3(f6.3,a),i5,a,i5)')ch10,&
325 &        ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), bands ',&
326 &        efmasdeg(ikpt)%degs_bounds(1,ideg),' through ',efmasdeg(ikpt)%degs_bounds(2,ideg)
327      endif
328 
329      sign_warn=.false.
330      do iband=1,deg_dim
331        if(saddle_warn(iband)) then
332          write(ab_out,'(a,i5,a)') ' Band ',efmasdeg(ikpt)%degs_bounds(1,ideg)+iband-1,&
333 &          ' SADDLE POINT - Frohlich effective mass and ZPR cannot be defined. '
334          sign_warn=.true.
335        else
336          m_avg_frohlich(iband) = DSIGN(m_avg_frohlich(iband),m_avg(iband))
337          zpr_frohlich_avg(iband) = -DSIGN(zpr_frohlich_avg(iband),m_avg(iband))
338          write(ab_out,'(a,i5,a,f14.10)') &
339 &          ' Band ',efmasdeg(ikpt)%degs_bounds(1,ideg)+iband-1,&
340 &          ' Angular average effective mass for Frohlich model (<m**0.5>)**2= ',m_avg_frohlich(iband)
341        endif
342        if(start_eigf3d_pos(iband) .neqv. start_eigf3d_pos(1))then
343          sign_warn=.true.
344        endif
345      enddo
346 
347      if(sign_warn .eqv. .false.)then
348        zpr_frohlich = four*pi* two**(-half) * (sum(zpr_frohlich_avg(1:deg_dim))/deg_dim) / cryst%ucvol
349        write(ab_out,'(2a)')&
350 &       ' Angular and band average effective mass and ZPR for Frohlich model.'
351        write(ab_out,'(a,es16.6)') &
352 &       ' Value of     (<<m**0.5>>)**2 = ',(sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim)**2
353        write(ab_out,'(a,es16.6)') &
354 &       ' Absolute Value of <<m**0.5>> = ', sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim
355        write(ab_out,'(a,es16.6,a,es16.6,a)') &
356 &       ' ZPR from Frohlich model      = ',zpr_frohlich,' Ha=',zpr_frohlich*Ha_eV,' eV'
357      else
358        write(ab_out,'(a)')& 
359 &        ' Angular and band average effective mass for Frohlich model cannot be defined because of a sign problem.'
360      endif
361 
362      ABI_DEALLOCATE(eig2_diag_cart)
363      ABI_DEALLOCATE(f3d)
364      ABI_DEALLOCATE(m_avg)
365      ABI_DEALLOCATE(m_avg_frohlich)
366      ABI_DEALLOCATE(zpr_frohlich_avg)
367      ABI_DEALLOCATE(eigenval)
368      ABI_DEALLOCATE(saddle_warn)
369      ABI_DEALLOCATE(start_eigf3d_pos)
370 
371    enddo ! ideg
372  enddo ! ikpt
373 
374  ABI_DEALLOCATE(unit_qdir)
375  ABI_DEALLOCATE(weight_qdir)
376  ABI_DEALLOCATE(polarity_qdir)
377  ABI_DEALLOCATE(proj_polarity_qdir)
378  ABI_DEALLOCATE(phfrq_qdir)
379  ABI_DEALLOCATE(dielt_qdir)
380  ABI_DEALLOCATE(zpr_q0_phononfactor_qdir)
381  ABI_DEALLOCATE(frohlich_phononfactor_qdir)
382 
383  end subroutine frohlichmodel
384 
385 end module m_frohlichmodel