TABLE OF CONTENTS


ABINIT/m_thmeig [ Modules ]

[ Top ] [ Modules ]

NAME

  m_thmeig

FUNCTION

 Calculate thermal corrections to the eigenvalues.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (PB, XG, GA)
  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

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_thmeig
28 
29  use defs_basis
30  use m_abicore
31  use m_tetrahedron
32  use m_errors
33  use m_ddb
34  use m_ddb_hdr
35  use m_xmpi
36  use m_sort
37 
38  use m_geometry,       only : mkrdim, xred2xcart, metric
39  use m_symfind,        only : symfind, symlatt
40  use m_symtk,          only : mati3inv, matr3inv, symatm
41  use m_crystal,        only : crystal_t
42  use m_io_tools,       only : open_file
43  use m_dynmat,         only : asria_corr, dfpt_phfrq
44  use m_anaddb_dataset, only : anaddb_dataset_type
45  use m_pawtab,         only : pawtab_type,pawtab_nullify,pawtab_free
46  use m_kpts,           only : getkgrid
47 
48  implicit none
49 
50  private

m_thmeig/outg2f [ Functions ]

[ Top ] [ m_thmeig ] [ Functions ]

NAME

 outg2f

FUNCTION

  Output g2f function to file. FIXME: Paul, please explain what g2f is.
  Probably a variant on the Eliashberg spectral function a2F

INPUTS

OUTPUT

  only write

PARENTS

      thmeig

CHILDREN

SOURCE

1152 subroutine outg2f(deltaene,enemin,enemax,filnam,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f)
1153 
1154 
1155 !This section has been created automatically by the script Abilint (TD).
1156 !Do not modify the following lines by hand.
1157 #undef ABI_FUNC
1158 #define ABI_FUNC 'outg2f'
1159 !End of the abilint section
1160 
1161  implicit none
1162 
1163 !Arguments ------------------------------------
1164 !scalars
1165  integer,intent(in) :: mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f
1166  character(len=fnlen),intent(in) :: filnam
1167  real(dp) :: deltaene,enemin,enemax,g2fsmear
1168 !arrays
1169  real(dp) :: g2f(mband,nkpt,nene),kpnt(3,nkpt,nqpt)
1170 
1171 !Local variables-------------------------------
1172 !scalars
1173  integer :: iband,ikpt,iomega,iost
1174  real(dp) :: omega
1175  character(len=fnlen) :: outfile
1176  character(len=500) :: message
1177 !arrays
1178 
1179 ! *************************************************************************
1180 
1181 !output the g2f
1182    outfile = trim(filnam) // '_G2F'
1183    open (unit=unit_g2f,file=outfile,status='unknown',iostat=iost)
1184    if (iost /= 0) then
1185      write (message,'(3a)')' thmeig : ERROR- opening file ',trim(outfile),' as new'
1186      MSG_ERROR(message)
1187    end if
1188 
1189    write(std_out,*) ' g2f function'
1190    write (unit_g2f,'(a)') '#'
1191    write (unit_g2f,'(a)') '# ABINIT package : g2f file'
1192    write (unit_g2f,'(a)') '#'
1193    write (unit_g2f,'(a,I10)') '#     number of qpoints integrated over : ', nqpt
1194    write (unit_g2f,'(a,I10)') '#     number of energy points : ', nene
1195    write (unit_g2f,'(a,E16.6,a,E16.6,a)') '#       between omega_min = ', enemin, &
1196 &   ' Ha and omega_max = ', enemax, ' Ha'
1197    if(telphint==1)then
1198      write (unit_g2f,'(a,E16.6)') '#   and the smearing width for gaussians is ', g2fsmear
1199      write (unit_g2f,'(a)') '#'
1200    end if
1201    if(telphint==0)then
1202      write (unit_g2f,'(a,I10)') '#   number of tetrahedrons', ntetra
1203      write (unit_g2f,'(a)') '#'
1204    end if
1205 
1206 !Write only the a2f function for the first K point
1207 !ikpt=1
1208    do ikpt=1,nkpt
1209      write(unit_g2f,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
1210      do iband=1,mband
1211        write(unit_g2f,*) 'band :', iband
1212        omega = enemin
1213        do iomega=1,nene
1214          write (unit_g2f,*) omega*Ha_eV*1000, g2f(iband, ikpt,iomega)
1215          omega=omega+deltaene
1216        end do
1217      end do
1218    end do
1219 
1220    close (unit=unit_g2f)
1221 
1222 end subroutine outg2f

m_thmeig/outphdos [ Functions ]

[ Top ] [ m_thmeig ] [ Functions ]

NAME

 outphdos

FUNCTION

  Print out phonon density of states

INPUTS

   deltaene = step on energy/frequency grid, in Hartree
   dos_phon = phonon DOS calculated on a grid
   enemin = minimal frequency
   enemax = maximal frequency
   filnam = file name for output to disk
   g2fsmear = smearing width
   nene = number of points on energy axis
   nqpt = number of q-points
   ntetra = number of tetrahedra, if tetrahedron interpolation is used
   telphint = flag for el-phonon interpolation method (to indicate Gaussian or tetrahedron integration)
   unit_phdos = unit for phonon DOS output

OUTPUT

  only write

SIDE EFFECTS

NOTES

   FIXME
   overcomplete inputs. Eliminate unit_phdos (just filnam) and deltaene (gotten from max-min/nene)

PARENTS

      thmeig

CHILDREN

SOURCE

1051 subroutine outphdos(deltaene,dos_phon,enemin,enemax,filnam,g2fsmear,nene,nqpt,ntetra,telphint,unit_phdos)
1052 
1053 
1054 !This section has been created automatically by the script Abilint (TD).
1055 !Do not modify the following lines by hand.
1056 #undef ABI_FUNC
1057 #define ABI_FUNC 'outphdos'
1058 !End of the abilint section
1059 
1060  implicit none
1061 
1062 !Arguments ------------------------------------
1063 !scalars
1064  integer,intent(in) :: nene,nqpt,ntetra,telphint,unit_phdos
1065  character(len=fnlen),intent(in) :: filnam
1066  real(dp) :: deltaene,enemin,enemax,g2fsmear
1067 !arrays
1068  real(dp) :: dos_phon(nene)
1069 
1070 !Local variables-------------------------------
1071 !scalars
1072  integer :: iomega,iost,step10
1073  real(dp) :: dos_effective,omega
1074  character(len=fnlen) :: outfile
1075  character(len=500) :: message
1076 !arrays
1077 
1078 ! *************************************************************************
1079 
1080    outfile = trim(filnam) // '_PDS'
1081    write(message, '(3a)')ch10,&
1082 &   ' Will write phonon DOS in file ',trim(outfile)
1083    call wrtout(ab_out,message,'COLL')
1084    call wrtout(std_out,message,'COLL')
1085 
1086    write(message, '(4a)')ch10,&
1087 &   ' For checking purposes, write ten values in the present file.',ch10,&
1088 &   '       Index    Energy (in Ha)      DOS '
1089    call wrtout(ab_out,message,'COLL')
1090    call wrtout(std_out,message,'COLL')
1091 
1092    open (unit=unit_phdos,file=outfile,status='replace',iostat=iost)
1093    if (iost /= 0) then
1094      write (message,'(3a)')' Opening file ',trim(outfile),' as new'
1095      MSG_ERROR(message)
1096    end if
1097 
1098    write (unit_phdos,'(a)') '#'
1099    write (unit_phdos,'(a)') '# ABINIT package : phonon DOS file'
1100    write (unit_phdos,'(a)') '#'
1101    write (unit_phdos,'(a,i10)') '#   Number of Qpoints integrated over : ', nqpt
1102    write (unit_phdos,'(a,i10)') '#   Number of energy points : ', nene
1103    write (unit_phdos,'(a,es16.6,a,es16.6,a)') '#   between omega_min = ', enemin, &
1104 &   ' Ha and omega_max = ', enemax, ' Ha'
1105    if(telphint==1)then
1106      write (unit_phdos,'(a,es16.6)') '#   The smearing width for gaussians is ', g2fsmear
1107    end if
1108    if(telphint==0)then
1109      write (unit_phdos,'(a,i10)') '#   Number of tetrahedrons', ntetra
1110    end if
1111    write (unit_phdos,'(a)') '#'
1112    write (unit_phdos,'(a)') '#      Index    Energy (in Ha)      DOS '
1113 
1114    omega = enemin
1115    do iomega=1,nene
1116      dos_effective=dos_phon(iomega)
1117      if(abs(dos_effective)<tol16)then
1118        dos_effective=zero
1119      end if
1120      step10=nene/10
1121      if(mod(iomega,step10)==1)write (std_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective
1122      if(mod(iomega,step10)==1)write (ab_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective
1123      write (unit_phdos, '(i10,es18.6,es18.6)')iomega, omega, dos_effective
1124      omega=omega+deltaene
1125    end do
1126 
1127    close (unit=unit_phdos)
1128 
1129  end subroutine outphdos

m_thmeig/thmeig [ Functions ]

[ Top ] [ m_thmeig ] [ Functions ]

NAME

 thmeig

FUNCTION

 This routine calculates the thermal corrections to the eigenvalues.
 The output is this quantity for the input k point.

INPUTS

  elph_base_name = root filename for outputs
  eig2_filnam = name of the eig2 database file
  comm=MPI communicator

OUTPUT

PARENTS

      anaddb

CHILDREN

SOURCE

  81 subroutine thmeig(inp, ddb, crystal, &
  82 &                 elph_base_name, eig2_filnam, ddbun, iout, &
  83 &                 natom, mpert, msize, d2asr, &
  84 &                 comm)
  85 
  86 
  87 !This section has been created automatically by the script Abilint (TD).
  88 !Do not modify the following lines by hand.
  89 #undef ABI_FUNC
  90 #define ABI_FUNC 'thmeig'
  91 !End of the abilint section
  92 
  93  implicit none
  94 
  95 !Arguments ------------------------------------
  96 !scalars
  97  integer,intent(inout) :: natom
  98  integer,intent(in) :: mpert,msize
  99  integer,intent(in) :: comm
 100  character(len=*),intent(in) :: elph_base_name, eig2_filnam
 101  integer,intent(in) :: ddbun,iout
 102  type(crystal_t), intent(inout) :: crystal
 103  type(anaddb_dataset_type),intent(inout) :: inp
 104  type(ddb_type),intent(inout) :: ddb
 105 !arrays
 106  real(dp),intent(inout) :: d2asr(2,3,natom,3,natom)
 107 
 108 
 109 !Local variables-------------------------------
 110 !scalars
 111  integer,parameter :: msppol=2,master=0,bcorr0=0
 112  integer :: msym
 113  integer :: nkpt,mband,ntypat
 114  integer :: usepaw,natifc
 115  integer :: nsym,occopt,nblok2
 116  integer :: ntemper,telphint,thmflag
 117  integer :: brav,chksymbreak,found,gqpt,iatom1,iatom2,iband,iblok,iblok2,idir1,idir2,ii,ikpt,ilatt,imod,index
 118  integer :: iomega,iqpt,iqpt1,iqpt2,iqpt2_previous,iqpt3,iscf_fake,itemper
 119  integer :: mpert_eig2,msize2,nene,ng2f,nqshft,nsym_new,unit_g2f,nqpt,nqpt_computed,qptopt,rftyp
 120 !integer :: mqpt,nqpt2,option
 121  integer :: unit_phdos,unitout
 122  integer :: isym
 123  integer :: nptsym,use_inversion
 124  integer :: ierr
 125  real(dp) :: ucvol
 126  real(dp) :: g2fsmear,temperinc,tempermin
 127  real(dp) :: bosein,deltaene,det,domega,enemax,enemin,fact2i,fact2r,factr
 128  real(dp) :: gaussfactor,gaussprefactor,gaussval,invdet,omega,omega_max,omega_min,qnrm,qptrlen
 129  real(dp) :: rcvol,tmp,tol,vec1i,vec1r,vec2i,vec2r,veci,vecr,xx
 130  real(dp) :: tolsym,tolsym8  !new
 131  character(len=500) :: message
 132  character(len=fnlen) :: outfile
 133  type(ddb_type) :: ddb_eig2
 134  type(ddb_hdr_type) :: ddb_hdr
 135 !arrays
 136  ! FIXME now these must be allocated
 137  integer :: ngqpt(9),qptrlatt(3,3),rfelfd(4),rfphon(4),rfstrs(4),vacuum(3)
 138  integer :: bravais(11)
 139  integer,allocatable :: typat(:),atifc(:)
 140  integer,allocatable :: symrel(:,:,:),symrec(:,:,:)
 141  integer,allocatable :: indsym(:,:,:)
 142  integer,allocatable :: indqpt(:)
 143  integer,allocatable :: symafm(:),symafm_new(:)
 144  integer,allocatable :: carflg_eig2(:,:,:,:)
 145  integer,allocatable :: ptsymrel(:,:,:),symrel_new(:,:,:)
 146 !integer,allocatable :: symrec_new(:,:,:)
 147  real(dp) :: rprim(3,3),gprim(3,3),rmet(3,3),gmet(3,3)
 148  real(dp) :: acell(3)
 149  real(dp) :: diff_qpt(3)
 150  real(dp) :: gprimd(3,3),mesh(3,3)
 151  real(dp) :: qlatt(3,3),qphnrm(3),qpt_search(3,3)
 152  real(dp) :: rprimd(3,3),shiftq(3,210),tempqlatt(3)
 153  real(dp) :: dummy(0),dummy2(0,0)
 154  real(dp),allocatable :: xcart(:,:),xred(:,:)
 155  real(dp),allocatable :: amu(:),zion(:)
 156  real(dp),allocatable :: tnons(:,:)
 157  real(dp),allocatable :: deigi(:,:), deigr(:,:), multi(:,:), multr(:,:)
 158  real(dp),allocatable :: dwtermi(:,:), dwtermr(:,:)
 159  real(dp),allocatable :: slope(:,:,:),thmeigen(:,:,:),zeropoint(:,:,:)
 160  real(dp),allocatable :: displ(:)
 161  real(dp),allocatable :: dos_phon(:),dtweightde(:,:),d2cart(:,:)
 162  real(dp),allocatable :: eigvec(:,:,:,:),eigval(:,:),g2f(:,:,:),intweight(:,:,:)
 163  real(dp),allocatable :: indtweightde(:,:,:),tmpg2f(:,:,:),tmpphondos(:),total_dos(:),tweight(:,:)
 164  real(dp),allocatable :: phfreq(:,:)
 165  real(dp),allocatable :: blkval2(:,:,:,:),blkval2gqpt(:,:,:,:),kpnt(:,:,:)
 166  real(dp),allocatable :: dedni(:,:,:,:),dednr(:,:,:,:)
 167  real(dp),allocatable :: eigen_in(:)
 168  real(dp),allocatable :: qpt_full(:,:),qptnrm(:)
 169  real(dp),allocatable :: spqpt(:,:),tnons_new(:,:),spinat(:,:)
 170  real(dp),allocatable :: wghtq(:)
 171 
 172  type(t_tetrahedron) :: tetrahedra
 173  character(len=80) :: errstr
 174 
 175 ! *********************************************************************
 176 
 177 !DEBUG
 178 ! write(std_out,*)'-thmeig : enter '
 179 !call flush(6)
 180 !ENDDEBUG
 181 
 182  ! Only master works for the time being
 183  if (xmpi_comm_rank(comm) /= master) return
 184 
 185  write(message,'(83a)') ch10,('=',ii=1,80),ch10,&
 186 & ' Computation of the electron-phonon changes to the electronic eigenenergies '
 187  call wrtout(ab_out,message,'COLL')
 188  call wrtout(std_out,message,'COLL')
 189 
 190 
 191 !=========================================================================
 192 !0) Initializations
 193 !=========================================================================
 194 
 195 
 196  g2fsmear = inp%a2fsmear
 197 
 198  telphint = inp%telphint
 199  temperinc = inp%temperinc
 200  tempermin = inp%tempermin
 201  thmflag = inp%thmflag
 202 
 203  ntemper = inp%ntemper
 204  natifc = inp%natifc
 205 
 206 !Open Derivative DataBase then r/w Derivative DataBase preliminary information.
 207 
 208  write(std_out, '(a)' )  '- thmeig: Initialize the second-order electron-phonon file with name :'
 209  write(std_out, '(a,a)' )'-         ',trim(eig2_filnam)
 210 
 211  call ddb_hdr_open_read(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION)
 212 
 213  mband = ddb_hdr%mband
 214  nkpt = ddb_hdr%nkpt
 215  ntypat = ddb_hdr%ntypat
 216 
 217  msym = ddb_hdr%msym
 218  nblok2 = ddb_hdr%nblok
 219  nsym = ddb_hdr%nsym
 220  occopt = ddb%occopt
 221  usepaw = ddb_hdr%usepaw
 222 
 223  ABI_ALLOCATE(typat, (natom))
 224  ABI_ALLOCATE(atifc, (natom))
 225  ABI_ALLOCATE(zion, (ntypat))
 226  ABI_ALLOCATE(amu, (ntypat))
 227 
 228  ABI_ALLOCATE(xcart,(3,natom))
 229  ABI_ALLOCATE(xred,(3,natom))
 230 
 231  ABI_ALLOCATE(symafm, (nsym))
 232  ABI_ALLOCATE(spinat,(3,natom))
 233 
 234  ABI_ALLOCATE(symrel, (3,3,nsym))
 235  ABI_ALLOCATE(symrec, (3,3,nsym))
 236  ABI_ALLOCATE(tnons, (3,nsym))
 237  ABI_ALLOCATE(indsym, (4,nsym,natom))
 238 
 239  ABI_ALLOCATE(deigi, (mband,nkpt))
 240  ABI_ALLOCATE(deigr, (mband,nkpt))
 241  ABI_ALLOCATE(dwtermi, (mband,nkpt))
 242  ABI_ALLOCATE(dwtermr, (mband,nkpt))
 243  ABI_ALLOCATE(multi, (mband,nkpt))
 244  ABI_ALLOCATE(multr, (mband,nkpt))
 245  ABI_ALLOCATE(slope, (2,mband,nkpt))
 246  ABI_ALLOCATE(thmeigen, (2,mband,nkpt))
 247  ABI_ALLOCATE(zeropoint, (2,mband,nkpt))
 248 
 249 !At present, only atom-type perturbations are allowed for eig2 type matrix elements.
 250  mpert_eig2=natom
 251  msize2=3*mpert_eig2*3*mpert_eig2
 252 
 253  call ddb_malloc(ddb_eig2,msize2,nblok2,natom,ntypat)
 254 
 255  ABI_ALLOCATE(blkval2,(2,msize2,mband,nkpt))
 256  ABI_ALLOCATE(blkval2gqpt,(2,msize2,mband,nkpt))
 257 
 258  ABI_ALLOCATE(eigvec,(2,3,natom,3*natom))
 259  ABI_ALLOCATE(phfreq,(3*natom,ddb%nblok))
 260 
 261  atifc = inp%atifc
 262 
 263  !amu = ddb%amu
 264  amu(:) = ddb_hdr%amu(1:ntypat)
 265  typat(:) = ddb_hdr%typat(1:natom)
 266  zion(:) = ddb_hdr%zion(1:ntypat)
 267  symrel(:,:,1:nsym) = ddb_hdr%symrel(:,:,1:nsym)
 268  tnons(:,1:nsym) = ddb_hdr%tnons(:,1:nsym)
 269 
 270  xred(:,:) = ddb_hdr%xred(:,:)
 271 
 272  symafm(:) = ddb_hdr%symafm(1:nsym)
 273  spinat(:,:) = ddb_hdr%spinat(:,1:natom)
 274 
 275  !symrel = ddb_hdr%symrel  ! out
 276  !tnons = ddb_hdr%tnons  ! out
 277 
 278  !acell = ddb%acell
 279  !natom = ddb_hdr%natom
 280  acell = ddb_hdr%acell
 281  rprim = ddb_hdr%rprim
 282 
 283 !Compute different matrices in real and reciprocal space, also
 284 !checks whether ucvol is positive.
 285  call mkrdim(acell,rprim,rprimd)
 286  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 287 
 288 !Obtain reciprocal space primitive transl g from inverse trans of r
 289  call matr3inv(rprim,gprim)
 290 
 291 !Generate atom positions in cartesian coordinates
 292  call xred2xcart(natom,rprimd,xcart,xred)
 293 
 294 !Transposed inversion of the symmetry matrices, for use in
 295 !the reciprocal space
 296  do isym=1,nsym
 297    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
 298  end do
 299 
 300 !SYMATM generates for all the atoms and all the symmetries, the atom
 301 !on which the referenced one is sent and also the translation bringing
 302 !back this atom to the referenced unit cell
 303  tolsym8=tol8
 304  call symatm(indsym,natom,nsym,symrec(:,:,1:nsym),tnons(:,1:nsym),tolsym8,typat,xred)
 305 
 306 !Check the correctness of some input parameters,
 307 !and perform small treatment if needed.
 308  call chkin9(atifc,natifc,natom)
 309 
 310  blkval2gqpt(:,:,:,:)=zero
 311 
 312  ABI_ALLOCATE(carflg_eig2,(3,mpert_eig2,3,mpert_eig2))
 313  ABI_ALLOCATE(kpnt,(3,nkpt,1))
 314 
 315  ! Copy a bunch of stuff back into crystal (to retain old behavior)
 316  ! TODO comment these: doesnt make a difference
 317  crystal%xcart = xcart
 318  crystal%ucvol = ucvol
 319  crystal%zion = zion
 320  crystal%gmet = gmet
 321  crystal%rmet = rmet
 322  crystal%nsym = nsym
 323  crystal%symrel = symrel
 324  crystal%symrec = symrec
 325  crystal%tnons = tnons
 326  crystal%indsym = indsym
 327 
 328 
 329 !=========================================================================
 330 !1) Take care of the Gamma point for thmflag=3, 5 or 7
 331 !=========================================================================
 332 
 333  if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
 334    found=0
 335    do iblok2=1,nblok2
 336 
 337      call read_blok8(ddb_eig2,iblok2,mband,mpert_eig2,msize2,&
 338 &     nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1))
 339 
 340 
 341      qnrm = ddb_eig2%qpt(1,iblok2)*ddb_eig2%qpt(1,iblok2)+ &
 342 &     ddb_eig2%qpt(2,iblok2)*ddb_eig2%qpt(2,iblok2)+ &
 343 &     ddb_eig2%qpt(3,iblok2)*ddb_eig2%qpt(3,iblok2)
 344      if(qnrm < DDB_QTOL) then
 345        blkval2gqpt(:,:,:,:) = blkval2(:,:,:,:)
 346        gqpt=iblok2
 347        write(std_out,*)'-thmeig: found Gamma point in EIG2 DDB, blok number ',iblok2
 348        found=1
 349        exit
 350      end if
 351    end do
 352 
 353    if(found==0)then
 354      write(message,'(a,i3,2a)')&
 355 &     'Was unable to find the blok for Gamma point in EIG2 DDB file, while thmflag= ',thmflag,ch10,&
 356 &     'Action: compute the contribution from Gamma, and merge it in your EIG2 DDB file.'
 357      MSG_ERROR(message)
 358    end if
 359 
 360 !  Put blkval2gqpt in cartesian coordinates
 361    call carttransf(ddb_eig2%flg,blkval2gqpt,carflg_eig2,gprimd,gqpt,mband,&
 362 &   mpert_eig2,msize2,natom,nblok2,nkpt,rprimd)
 363 
 364  end if
 365 
 366  close(ddbun)
 367 
 368 !=========================================================================
 369 !2) Calculation of dE(n,k)/dn(Q,j) : consider all q and modes
 370 !=========================================================================
 371 
 372  if(thmflag==3 .or. thmflag==4)then
 373 
 374 
 375 !  Use the first list of q wavevectors
 376    nqpt=inp%nph1l
 377    ABI_ALLOCATE(spqpt,(3,nqpt))
 378    do iqpt=1,inp%nph1l
 379      spqpt(:,iqpt)=inp%qph1l(:,iqpt)/inp%qnrml1(iqpt)
 380    end do
 381    ABI_ALLOCATE(wghtq,(nqpt))
 382    wghtq(:)=one/nqpt
 383 
 384  else if(thmflag>=5 .and. thmflag<=8)then
 385 
 386 !  Generates the q point grid
 387    ngqpt(1:3)=inp%ngqpt(1:3)
 388    nqshft=inp%nqshft
 389    qptrlatt(:,:)=0
 390    qptrlatt(1,1)=ngqpt(1)
 391    qptrlatt(2,2)=ngqpt(2)
 392    qptrlatt(3,3)=ngqpt(3)
 393 
 394    ABI_ALLOCATE(ptsymrel,(3,3,msym))
 395    ABI_ALLOCATE(symafm_new,(msym))
 396    ABI_ALLOCATE(symrel_new,(3,3,msym))
 397    ABI_ALLOCATE(tnons_new,(3,msym))
 398    if(thmflag==7 .or. thmflag==8) then
 399 !    Re-generate symmetry operations from the lattice and atomic coordinates
 400      tolsym=tol8
 401      call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
 402      use_inversion=1
 403      call symfind(0,(/zero,zero,zero/),gprimd,0,msym,natom,0,nptsym,nsym_new,0,0,&
 404 &     ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred)
 405      write(std_out,*)' thmeig : found ',nsym_new,' symmetries ',ch10
 406      qptopt=1
 407    else
 408      nsym_new=1
 409      symrel_new(:,:,1)=0 ; symrel_new(1,1,1)=1 ; symrel_new(2,2,1)=1 ; symrel_new(3,3,1)=1
 410      tnons_new(:,1)=zero
 411      symafm_new(1)=1
 412      qptopt=3
 413    end if
 414 
 415    brav=inp%brav
 416 
 417    if(abs(brav)/=1)then
 418      message = ' The possibility to have abs(brav)/=1 for thmeig was disabled.'
 419      MSG_ERROR(message)
 420    end if
 421 
 422 !  Prepare to compute the q-point grid in the ZB or IZB
 423    iscf_fake=5 ! Need the weights
 424    chksymbreak=0
 425    vacuum=0
 426    shiftq(:,1:nqshft)=inp%q1shft(:,1:nqshft)
 427 !  Compute the final number of q points
 428    call getkgrid(chksymbreak,0,iscf_fake,dummy2,qptopt,qptrlatt,qptrlen,&
 429 &   nsym_new,0,nqpt,nqshft,nsym_new,rprimd,&
 430 &   shiftq,symafm_new,symrel_new,vacuum,dummy)
 431    ABI_ALLOCATE(spqpt,(3,nqpt))
 432    ABI_ALLOCATE(wghtq,(nqpt))
 433    call getkgrid(chksymbreak,iout,iscf_fake,spqpt,qptopt,qptrlatt,qptrlen,&
 434 &   nsym_new,nqpt,nqpt_computed,nqshft,nsym_new,rprimd,&
 435 &   shiftq,symafm_new,symrel_new,vacuum,wghtq)
 436 
 437    ABI_DEALLOCATE(ptsymrel)
 438    ABI_DEALLOCATE(symafm_new)
 439    ABI_DEALLOCATE(symrel_new)
 440    ABI_DEALLOCATE(tnons_new)
 441 
 442  end if
 443 
 444  call ddb_hdr_free(ddb_hdr)
 445 
 446 
 447  write(message,'(a,a)')ch10,' thmeig : list of q wavevectors, with integration weights '
 448  call wrtout(ab_out,message,'COLL')
 449  call wrtout(std_out,message,'COLL')
 450  do iqpt=1,nqpt
 451    write(message,'(i6,3es16.6,es20.6)')iqpt,spqpt(:,iqpt),wghtq(iqpt)
 452    call wrtout(ab_out,message,'COLL')
 453    call wrtout(std_out,message,'COLL')
 454  end do
 455 
 456  if(.not.allocated(indqpt))then
 457    ABI_ALLOCATE(indqpt,(nqpt))
 458  end if
 459  ABI_ALLOCATE(dedni,(mband,nkpt,3*natom,nqpt))
 460  ABI_ALLOCATE(dednr,(mband,nkpt,3*natom,nqpt))
 461  ABI_ALLOCATE(eigen_in,(nqpt))
 462  ABI_ALLOCATE(qpt_full,(3,nqpt))
 463  ABI_ALLOCATE(qptnrm,(nqpt))
 464 
 465  dednr(:,:,:,:) = zero
 466  dedni(:,:,:,:) = zero
 467 
 468 !!Prepare the reading of the EIG2 files
 469  call ddb_hdr_open_read(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION, msym=msym)
 470 
 471  call ddb_hdr_free(ddb_hdr)
 472 
 473 !iqpt2 will be the index of the q point bloks inside the EIG2 file
 474  iqpt2=0
 475 
 476 !Sum on all phonon wavevectors and modes
 477  do iqpt=1,nqpt
 478 
 479 !  Finding the target wavevector in DDB file
 480    qpt_search(:,:)=0.0d0
 481    qpt_search(:,1)=spqpt(:,iqpt)
 482    qphnrm(:)=one
 483    rfphon(1:2)=1
 484 !  NOTE : at present, no LO-TO splitting included !!!
 485    rfelfd(1:2)=0
 486    rfstrs(1:2)=0
 487    rftyp=1
 488 
 489    write(std_out,'(a,3es16.6)' )' Looking for spqpt=',qpt_search(:,1)
 490 
 491    call gtblk9(ddb,iblok,qpt_search,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
 492 
 493    if(iblok==0) then
 494      write(message,'(a,3es16.6,2a)')&
 495 &     'Was unable to find in DDB file, the blok for point ',spqpt(:,iqpt),ch10,&
 496 &     'Action: compute the contribution from this point, and merge it in your DDB file.'
 497      MSG_ERROR(message)
 498    end if
 499 
 500    ABI_ALLOCATE(d2cart,(2,msize))
 501 !  Copy the dynamical matrix in d2cart
 502    d2cart(:,1:msize)=ddb%val(:,:,iblok)
 503 
 504 !  Eventually impose the acoustic sum rule based on previously calculated d2asr
 505    !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
 506    if (inp%asr==1 .or. inp%asr==2 .or. inp%asr==5) then
 507      call asria_corr(inp%asr,d2asr,d2cart,mpert,natom)
 508    end if
 509 
 510 !  Calculation of the eigenvectors and eigenvalues
 511 !  of the dynamical matrix
 512    ABI_ALLOCATE(displ,(2*3*natom*3*natom))
 513    ABI_ALLOCATE(eigval,(3,natom))
 514    call dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,&
 515 &   mpert,msym,natom,nsym,ntypat,phfreq(:,iqpt),qphnrm(1),spqpt(:,iqpt),rprimd,inp%symdynmat,&
 516 &   symrel,symafm,typat,ucvol)
 517    ABI_DEALLOCATE(displ)
 518    ABI_DEALLOCATE(eigval)
 519    ABI_DEALLOCATE(d2cart)
 520 
 521 
 522 !  Read the next bloks to find the next q point.
 523    found=0 ; iqpt2_previous=iqpt2
 524    do while (iqpt2<nblok2)
 525      iqpt2=iqpt2+1
 526      call read_blok8(ddb_eig2,iqpt2,mband,mpert_eig2,msize2,&
 527 &     nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1))
 528      !write (300,*) 'blkval2 _bis_ in thmeig'
 529      !write (300,*) blkval2
 530      diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt)
 531      if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then
 532        found=1
 533        exit
 534      end if
 535    end do
 536 
 537 !  Usually, the q points come in the right order. However, this is not always the case...
 538    if(found==0)then
 539 
 540 !    If the EIG2 database file has to be read again, close it, then search for the right q point,
 541 !    from the beginning of the file
 542      close(ddbun)
 543 
 544      call ddb_hdr_open_read(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION, msym=msym)
 545      call ddb_hdr_free(ddb_hdr)
 546 
 547 !    And examine again the EIG2 file. Still, not beyond the previously examined value.
 548      found=0
 549      do iqpt2=1,iqpt2_previous
 550        call read_blok8(ddb_eig2,iqpt2,mband,mpert_eig2,msize2,&
 551 &       nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1))
 552        diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt)
 553        if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then
 554          found=1
 555          exit
 556        end if
 557      end do
 558 
 559      if(found==0)then
 560        write(message,'(a,3es16.6,2a)')&
 561 &       'Was unable to find in EIG2 DDB file, the blok for point ',spqpt(:,iqpt),ch10,&
 562 &       'Action: compute the contribution from this point, and merge it in your EIG2 DDB file.'
 563        MSG_ERROR(message)
 564      end if
 565 
 566    end if
 567 
 568 !  Put blkval2 in cartesian coordinates
 569    call carttransf(ddb_eig2%flg,blkval2,carflg_eig2,gprimd,iqpt,mband,&
 570 &   mpert_eig2,msize2,natom,nblok2,nkpt,rprimd)
 571 
 572    do imod=1,3*natom
 573 
 574 !    Calculate the derivative
 575      deigr(:,:) = zero
 576      deigi(:,:) = zero
 577      dwtermr(:,:)=zero
 578      dwtermi(:,:)=zero
 579      index=0
 580      do iatom1=1,natom
 581        do idir1=1,3
 582          do iatom2=1,natom
 583 !          Compute factor for SE term
 584            if(phfreq(imod,iqpt)<tol6)then
 585              factr = zero
 586            else
 587              factr=one/sqrt(amu(typat(iatom1))*amu(typat(iatom2)))/phfreq(imod,iqpt)/amu_emass
 588            end if
 589 
 590            do idir2=1,3
 591              index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1)))
 592 
 593 !            Compute products of polarization vectors
 594              vecr = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)+&
 595 &             eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod)
 596              veci = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)-&
 597 &             eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod)
 598 
 599              vec1r = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)+&
 600 &             eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod)
 601              vec1i = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)-&
 602 &             eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod)
 603 
 604              vec2r = eigvec(1,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)+&
 605 &             eigvec(2,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod)
 606              vec2i = eigvec(2,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)-&
 607 &             eigvec(1,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod)
 608 
 609 !            Compute factor for DW term
 610              if(phfreq(imod,iqpt)<tol6)then
 611                fact2r = zero
 612                fact2i = zero
 613              else
 614                fact2r = -wghtq(iqpt)*(vec1r/amu(typat(iatom1)) + vec2r/amu(typat(iatom2)))/phfreq(imod,iqpt)/&
 615 &               amu_emass/2 !/norm(idir1)/norm(idir2)
 616                fact2i = -wghtq(iqpt)*(vec1i/amu(typat(iatom1)) + vec2i/amu(typat(iatom2)))/phfreq(imod,iqpt)/&
 617 &               amu_emass/2 !/norm(idir1)/norm(idir2)
 618              end if
 619 
 620              multr(:,:) =(blkval2(1,index,:,:)*vecr - blkval2(2,index,:,:)*veci) !/(norm(idir1)*norm(idir2))
 621              multi(:,:) =(blkval2(1,index,:,:)*veci + blkval2(2,index,:,:)*vecr) !/(norm(idir1)*norm(idir2))
 622 
 623 
 624 !            Debye-Waller Term
 625              if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
 626                dwtermr(1:mband,1:nkpt)=dwtermr(1:mband,1:nkpt)+fact2r*blkval2gqpt(1,index,:,:)-fact2i*blkval2gqpt(2,index,:,:)
 627                dwtermi(1:mband,1:nkpt)=dwtermi(1:mband,1:nkpt)+fact2r*blkval2gqpt(2,index,:,:)+fact2i*blkval2gqpt(1,index,:,:)
 628              end if
 629 
 630 !            Self-energy Term (Fan)
 631              deigr(1:mband,1:nkpt) = deigr(1:mband,1:nkpt) + wghtq(iqpt)*factr*multr(1:mband,1:nkpt)
 632              deigi(1:mband,1:nkpt) = deigi(1:mband,1:nkpt) + wghtq(iqpt)*factr*multi(1:mband,1:nkpt)
 633 
 634            end do !idir2
 635          end do !iatom2
 636        end do !idir1
 637      end do !iatom1
 638 !    Eigenvalue derivative or broadening
 639      if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then
 640        dednr(1:mband,1:nkpt,imod,iqpt) = deigr(1:mband,1:nkpt) + dwtermr(1:mband,1:nkpt)
 641        dedni(1:mband,1:nkpt,imod,iqpt) = deigi(1:mband,1:nkpt) + dwtermi(1:mband,1:nkpt)
 642      else if(thmflag==4 .or. thmflag==6 .or. thmflag==8) then
 643        dednr(1:mband,1:nkpt,imod,iqpt) = pi*deigr(1:mband,1:nkpt)
 644        dedni(1:mband,1:nkpt,imod,iqpt) = pi*deigi(1:mband,1:nkpt)
 645      end if
 646 
 647    end do ! imod
 648  end do !iqpt
 649 
 650  close(ddbun)
 651 
 652 
 653 !=============================================================================
 654 !3) Evaluation of the Eliashberg type spectral function
 655 !and phonon DOS via gaussian broadning
 656 !=============================================================================
 657 
 658  if(telphint==1)then
 659    ng2f = 500  ! number of frequencies
 660    omega_min=zero
 661    omega_max=zero
 662    do iqpt=1,nqpt
 663      do imod=1,3*natom
 664        omega_min = min(omega_min,phfreq(imod,iqpt))
 665        omega_max = max(omega_max,phfreq(imod,iqpt))
 666      end do
 667    end do
 668 
 669    ABI_ALLOCATE(dos_phon,(ng2f))
 670    ABI_ALLOCATE(g2f,(mband,nkpt,ng2f))
 671    ABI_ALLOCATE(tmpg2f,(mband,nkpt,ng2f))
 672    ABI_ALLOCATE(tmpphondos,(ng2f))
 673 
 674    write(std_out,'(a,es13.6)') 'omega_min :', omega_min
 675    write(std_out,'(a,es13.6)') 'omega_max :', omega_max
 676    write(std_out,'(a,i8)') 'ng2f :', ng2f
 677 
 678    omega_max = omega_max + 0.1 * omega_max
 679    domega = (omega_max-omega_min)/(ng2f-one)
 680 
 681    gaussprefactor = sqrt(piinv) / g2fsmear
 682    gaussfactor = one / g2fsmear
 683 
 684    g2f(:,:,:) = zero
 685    dos_phon(:) = zero
 686 
 687    do iqpt=1,nqpt
 688      do imod=1,3*natom
 689        omega = omega_min
 690        tmpg2f(:,:,:) = zero
 691        tmpphondos(:) = zero
 692        do iomega=1,ng2f
 693          xx = (omega-phfreq(imod,iqpt))*gaussfactor
 694          gaussval = gaussprefactor*exp(-xx*xx)
 695          tmpg2f(:,:,iomega) = tmpg2f(:,:,iomega) + gaussval*dednr(:,:,imod,iqpt)
 696          tmpphondos(iomega) = tmpphondos(iomega) + gaussval
 697          omega = omega+domega
 698        end do
 699 
 700        g2f(:,:,:) = g2f(:,:,:) + tmpg2f(:,:,:)
 701        dos_phon(:) = dos_phon(:) + tmpphondos(:)
 702 
 703      end do !imod
 704    end do !iqpt
 705 
 706    dos_phon(:) = dos_phon(:) / nqpt
 707 
 708 !  output the g2f
 709    unit_g2f = 108
 710    call outg2f(domega,omega_min,omega_max,elph_base_name,g2f,g2fsmear,kpnt,mband,ng2f,nkpt,nqpt,1,telphint,unit_g2f)
 711 
 712 !  output the phonon DOS
 713    unit_phdos = 108
 714    call outphdos(domega,dos_phon,omega_min,omega_max,elph_base_name,g2fsmear,ng2f,nqpt,1,telphint,unit_g2f)
 715 
 716 
 717    ABI_DEALLOCATE(dos_phon)
 718    ABI_DEALLOCATE(g2f)
 719    ABI_DEALLOCATE(tmpg2f)
 720    ABI_DEALLOCATE(tmpphondos)
 721 
 722  end if !telphint
 723 
 724 !=======================================================================
 725 !4) Evaluation of the Eliashberg type spectral function
 726 !and phonon DOS via improved tetrahedron method
 727 !=======================================================================
 728 
 729  if(telphint==0)then
 730 
 731 !  make dimension-ful rprimd and gprimd for transformation of derivatives to cartesian coordinates.
 732    call mkrdim(acell,rprim,rprimd)
 733    call matr3inv(rprimd,gprimd)
 734 
 735 !  Q point Grid
 736    qpt_full(:,:) = ddb%qpt(1:3,:)
 737 
 738 !  Trivial Q point index
 739    do iqpt=1,nqpt
 740      indqpt(iqpt)=iqpt
 741      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)
 742    end do
 743 
 744 !  Build qlatt from scratch (for 5.7)
 745    tol = 0.1_dp
 746    ilatt = 0
 747    call sort_dp(nqpt,qptnrm,indqpt,tol)
 748 
 749    do iqpt1=1,nqpt-2
 750      mesh(1:3,1) = qpt_full(1:3,indqpt(iqpt1))
 751      do iqpt2=iqpt1+1,nqpt-1
 752        mesh(1:3,2)= qpt_full(1:3,indqpt(iqpt2))
 753        do iqpt3=iqpt2+1,nqpt
 754          mesh(1:3,3)= qpt_full(1:3,indqpt(iqpt3))
 755          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) &
 756 &         -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)
 757          invdet = one/det
 758          if (abs(nint(invdet))==nqpt .and. abs(invdet)-nqpt < tol) then
 759            ilatt = 1
 760            qlatt(:,:) = mesh(:,:)
 761            exit
 762          end if
 763        end do
 764        if(ilatt==1) exit
 765      end do
 766      if(ilatt==1) exit
 767    end do
 768 
 769 !  error message if qlatt not found and stop
 770    if(ilatt==0) then
 771      write(message, '(a,a)' ) &
 772 &     ' Could not find homogeneous basis vectors for Q point grid ',ch10
 773      call wrtout(std_out,message,'COLL')
 774      call wrtout(ab_out,message,'COLL')
 775      MSG_ERROR("Aborting now")
 776    end if
 777 
 778 !  test if qlatt is righthanded and possibly fixe it
 779    if(invdet < 0) then
 780      tempqlatt(:) = qlatt(:,2)
 781      qlatt(:,2) = qlatt(:,1)
 782      qlatt(:,1) = tempqlatt(:)
 783    end if
 784 
 785    write(std_out,*) 'qlatt',qlatt
 786 
 787 !  test if qlatt generates all Q points  TO DO
 788 
 789 
 790 
 791 !  Get tetrahedra, ie indexes of the full kpoints at their summits
 792    call init_tetra(indqpt,gprimd,qlatt,qpt_full,nqpt,&
 793 &   tetrahedra, ierr, errstr)
 794    ABI_CHECK(ierr==0,errstr)
 795 
 796    rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
 797 &   -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
 798 &   +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
 799 
 800 !  Calculate weights for phonon DOS
 801 !  Special precautions must be taking for Gamma point
 802 !  because of non-analytic term.
 803 !  Non-analyticity must be taken out and treated separatly.
 804 
 805    nene = 100     !nene=number of energies for DOS
 806    enemin = minval(phfreq)
 807    enemax = maxval(phfreq)
 808    deltaene = (enemax-enemin)/dble(nene-1)
 809 !  redefine enemin enemax to be at rounded multiples of deltaene
 810 !  enemin = elph_ds%fermie - dble(ifermi)*deltaene
 811 !  enemax = elph_ds%fermie + dble(nene-ifermi-1)*deltaene
 812 
 813    ABI_ALLOCATE(tweight,(nqpt,nene))
 814    ABI_ALLOCATE(dtweightde,(nqpt,nene))
 815    ABI_ALLOCATE(intweight,(3*natom,nqpt,nene))
 816    ABI_ALLOCATE(indtweightde,(3*natom,nqpt,nene))
 817 
 818    do iband=1,3*natom
 819      eigen_in(:) = phfreq(iband,:)
 820 
 821 !    calculate general integration weights at each irred kpoint 
 822 !    as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
 823      call get_tetra_weight(eigen_in,enemin,enemax,&
 824 &     one,nene,nqpt,tetrahedra,bcorr0,&
 825 &     tweight,dtweightde,xmpi_comm_self)
 826 
 827      intweight(iband,:,:) = tweight(:,:)
 828      indtweightde(iband,:,:) = dtweightde(:,:)
 829 
 830    end do !iband
 831 
 832 !  intdtweightse(nband,nqpt,nene) represents the weight in each energy bin for every kpt and every band
 833 !  So phonon DOS is calculated (neglecting the non-analyticity contribution for now !!!)
 834 
 835    ABI_ALLOCATE(total_dos,(nene))
 836    ABI_ALLOCATE(g2f,(mband,nkpt,nene))
 837 
 838    total_dos(:) = zero
 839    do iband=1,3*natom
 840      do iqpt=1,nqpt
 841        total_dos(:) = total_dos + indtweightde(iband,iqpt,:)
 842      end do
 843    end do
 844 
 845 !  For the g2f function
 846 !  Right now for one electronic band and one K point: dednr(1:mband,1:nkpt,imod,iqpt)
 847 !  Once again must pay close attention to the Gamma point
 848    g2f(:,:,:) = zero
 849    do ii=1,mband
 850      do ikpt=1,nkpt
 851        do iband=1,3*natom
 852          do iqpt=1,nqpt
 853            g2f(ii,ikpt,:) = g2f(ii,ikpt,:) + dednr(ii,ikpt,iband,iqpt) * indtweightde(iband,iqpt,:)
 854          end do
 855        end do
 856      end do
 857    end do
 858 
 859 !  output the g2f
 860    unit_g2f = 108
 861    call outg2f(deltaene,enemin,enemax,elph_base_name,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,tetrahedra%ntetra,telphint,unit_g2f)
 862 
 863 !  output the phonon DOS
 864    unit_phdos = 108
 865    call outphdos(deltaene,total_dos,enemin,enemax,elph_base_name,g2fsmear,nene,nqpt,tetrahedra%ntetra,telphint,unit_g2f)
 866 
 867    ABI_DEALLOCATE(tweight)
 868    ABI_DEALLOCATE(dtweightde)
 869    ABI_DEALLOCATE(intweight)
 870    ABI_DEALLOCATE(indtweightde)
 871    ABI_DEALLOCATE(total_dos)
 872    ABI_DEALLOCATE(g2f)
 873  end if !telphint
 874 
 875 !=======================================================================
 876 !5) direct evaluation of thermal corrections
 877 !=======================================================================
 878 
 879 !open TBS file
 880  outfile = trim(elph_base_name)//"_TBS"
 881  if (open_file(outfile,message,newunit=unitout,form='formatted',status='unknown') /= 0) then
 882    MSG_ERROR(message)
 883  end if
 884  write(unitout,'(a)')'thmeig: Thermal Eigenvalue corrections (eV)'
 885 
 886  slope(:,:,:) = zero
 887  zeropoint(:,:,:) = zero
 888 !Loop on temperatures
 889  do itemper= 1, ntemper
 890    tmp=tempermin+temperinc*float(itemper-1)
 891    thmeigen(:,:,:) = zero
 892 
 893 !  Sum on all phonon wavevectors and modes
 894    do iqpt=1,nqpt
 895      do imod=1,3*natom
 896 
 897 !      Bose-Einstein distribution
 898 ! jmb overflow with exp(). So, select bosein to be still significant wrt half
 899        if(phfreq(imod,iqpt)<tol6 .or. (phfreq(imod,iqpt)/(kb_HaK*tmp)) > -log(tol16))then
 900          bosein = zero
 901        else
 902          bosein = one/(exp(phfreq(imod,iqpt)/(kb_HaK*tmp))-one)
 903        end if
 904 
 905 !      Calculate total
 906        thmeigen(1,1:mband,1:nkpt) = thmeigen(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(bosein+half)
 907        thmeigen(2,1:mband,1:nkpt) = thmeigen(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(bosein+half)
 908 
 909        if(itemper==1)then
 910 !        Calculate slope of linear regime
 911          if(phfreq(imod,iqpt)<tol6)then
 912            slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt)
 913            slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt)
 914          else
 915            slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt))
 916            slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt))
 917          end if
 918 !        Calculate zero-point renormalization
 919          zeropoint(1,1:mband,1:nkpt) = zeropoint(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*half
 920          zeropoint(2,1:mband,1:nkpt) = zeropoint(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*half
 921 
 922        end if
 923      end do ! imod
 924    end do !iqpt
 925 
 926 !  Write temperature independent results
 927    if(itemper==1)then
 928      write(unitout,'(a)')'Temperature independent results (zero-point renormalization and slope)'
 929      do ikpt=1,nkpt
 930        write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
 931        do iband=1,mband
 932          write(unitout,'(4d22.14)') Ha_eV*zeropoint(1,iband,ikpt),Ha_eV*zeropoint(2,iband,ikpt),&
 933 &         Ha_eV*slope(1,iband,ikpt),Ha_eV*slope(2,iband,ikpt)
 934        end do
 935      end do
 936      write(unitout,'(a)')'Temperature dependent corrections'
 937    end if
 938 !  Write result in file for each temperature
 939    write(unitout,'(a,es10.3,a)')'T :', tmp,' K'
 940    do ikpt=1,nkpt
 941      write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1)
 942      do iband=1,mband
 943        write(unitout,'(2d22.14)') Ha_eV*thmeigen(1,iband,ikpt), Ha_eV*thmeigen(2,iband,ikpt)
 944      end do
 945    end do
 946  end do !itemper
 947 
 948  close(unitout)
 949 
 950 !Write temperature-independent results to the main output file
 951  write(iout,'(a)')' '
 952  write(iout,'(80a)') ('-',ii=1,80)
 953  write(iout,'(a)')' '
 954  write(iout,'(a)')' Electron-phonon change of electronic structure.'
 955  write(iout,'(a)')' The temperature-dependent values are written in the _TBS file.'
 956  write(iout,'(a)')' Here follows, for each electronic wavevector and band :'
 957  write(iout,'(a)')'      zero-point renormalisation (Ha) and linear slope (Ha/Kelvin)'
 958  do ikpt=1,nkpt
 959    write(iout,'(2a,i6,a,3es16.6)')ch10,' Kpt number ',ikpt,', with reduced coordinates :',kpnt(:,ikpt,1)
 960    do iband=1,mband
 961      write(iout,'(i6,2es20.6)') iband,zeropoint(1,iband,ikpt),slope(1,iband,ikpt)
 962    end do
 963  end do
 964 
 965  ABI_DEALLOCATE(typat)
 966  ABI_DEALLOCATE(atifc)
 967  ABI_DEALLOCATE(zion)
 968  ABI_DEALLOCATE(amu)
 969  ABI_DEALLOCATE(xcart)
 970  ABI_DEALLOCATE(xred)
 971  ABI_DEALLOCATE(symafm)
 972  ABI_DEALLOCATE(spinat)
 973  ABI_DEALLOCATE(symrel)
 974  ABI_DEALLOCATE(symrec)
 975  ABI_DEALLOCATE(indsym)
 976  ABI_DEALLOCATE(tnons)
 977  ABI_DEALLOCATE(deigi)
 978  ABI_DEALLOCATE(deigr)
 979  ABI_DEALLOCATE(dwtermi)
 980  ABI_DEALLOCATE(dwtermr)
 981  ABI_DEALLOCATE(multi)
 982  ABI_DEALLOCATE(multr)
 983  ABI_DEALLOCATE(slope)
 984  ABI_DEALLOCATE(thmeigen)
 985  ABI_DEALLOCATE(zeropoint)
 986 
 987  ABI_DEALLOCATE(dedni)
 988  ABI_DEALLOCATE(dednr)
 989  if(allocated(indqpt)) then
 990    ABI_DEALLOCATE(indqpt)
 991  end if
 992  ABI_DEALLOCATE(eigen_in)
 993  ABI_DEALLOCATE(qpt_full)
 994  ABI_DEALLOCATE(qptnrm)
 995  ABI_DEALLOCATE(wghtq)
 996  ABI_DEALLOCATE(spqpt)
 997  ABI_DEALLOCATE(eigvec)
 998  ABI_DEALLOCATE(phfreq)
 999 
1000  ABI_DEALLOCATE(blkval2)
1001  ABI_DEALLOCATE(blkval2gqpt)
1002  ABI_DEALLOCATE(kpnt)
1003  ABI_DEALLOCATE(carflg_eig2)
1004 
1005 
1006 
1007  call ddb_free(ddb_eig2)
1008 
1009  call destroy_tetra(tetrahedra)
1010 
1011 end subroutine thmeig