TABLE OF CONTENTS


ABINIT/m_harmonic_thermo [ Modules ]

[ Top ] [ Modules ]

NAME

 m_harmonic_thermo

FUNCTION

 This routine to calculate phonon density of states,
 thermodynamical properties, Debye-Waller factor, and atomic mean square velocity

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (CL, 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_harmonic_thermo
29 
30  use defs_basis
31  use m_errors
32  use m_abicore
33  use m_sortph
34  use m_xmpi
35 
36  use m_io_tools,       only : open_file
37  use m_symtk,          only : matr3inv
38  use m_dynmat,         only : gtdyn9
39  use m_geometry,       only : mkrdim
40  use m_crystal,        only : crystal_t
41  use m_anaddb_dataset, only : anaddb_dataset_type
42  use m_ifc,            only : ifc_type, ifc_fourq
43  use m_kpts,           only : smpbz
44  use m_symkpt,     only : symkpt
45 
46  implicit none
47 
48  private

m_harmonic_thermo/harmonic_thermo [ Functions ]

[ Top ] [ m_harmonic_thermo ] [ Functions ]

NAME

 harmonic_thermo

FUNCTION

 This routine to calculate phonon density of states,
 thermodynamical properties, Debye-Waller factor, and atomic mean square velocity

INPUTS

 Crystal<crystal_t>=data type gathering info on the crystalline structure.
 Ifc<ifc_type>=Object containing the interatomic force constants.
    %atmfrc(2,3,natom,3,natom,nrpt) = Interatomic Forces in real space
    %dyewq0(3,3,natom)=Ewald part of the dynamical matrix, at q=0.
    %rpt(3,nrpt)=canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!)
    %nrpt=number of R points in the Big Box
    %trans(3,natom)=atomic translations : xred = rcan + trans
    %wghatm(natom,natom,nrpt)=weights associated to a pair of atoms and to a R vector
 amu(ntypat)=mass of the atoms (atomic mass unit)
 anaddb_dtset= (derived datatype) contains all the input variables
 iout =unit number for output
 natom=number of atoms in the unit cell
 outfilename_radix=radix of anaddb output file name: append _THERMO for thermodynamic quantities
 comm=MPI communicator

OUTPUT

NOTES

 dosinc=increment between the channels for the phonon DOS in cm-1

PARENTS

      anaddb

CHILDREN

      end_sortph,ifc_fourq,matr3inv,mkrdim,smpbz,sortph,symkpt,wrtout

SOURCE

  95 subroutine harmonic_thermo(Ifc,Crystal,amu,anaddb_dtset,iout,outfilename_radix,comm,thmflag)
  96 
  97 
  98 !This section has been created automatically by the script Abilint (TD).
  99 !Do not modify the following lines by hand.
 100 #undef ABI_FUNC
 101 #define ABI_FUNC 'harmonic_thermo'
 102 !End of the abilint section
 103 
 104  implicit none
 105 
 106 !Arguments -------------------------------
 107 !scalars
 108  integer,intent(in) :: iout,comm
 109  integer,intent(in),optional :: thmflag
 110  character(len=*),intent(in) :: outfilename_radix
 111  type(anaddb_dataset_type),intent(in) :: anaddb_dtset
 112  type(crystal_t),intent(in) :: Crystal
 113  type(ifc_type),intent(in) :: Ifc
 114 !arrays
 115  real(dp),intent(in) :: amu(Crystal%ntypat)
 116 
 117 !Local variables -------------------------
 118 !scalars
 119  integer,parameter :: master=0
 120  integer :: convth,facbrv,iatom,ichan,icomp,igrid,iii,iiii,ij,ind
 121  integer :: iqpt2,isym,itemper,iwchan,jjj,mqpt2,nchan,ngrids,natom
 122  integer :: nqpt2,nspqpt,ntemper,nwchan,option,timrev
 123  integer :: thermal_unit
 124  integer :: bij_unit
 125  integer :: vij_unit
 126  integer :: nomega, iomega
 127  real(dp) :: change,cothx,diffbb,dosinc,expm2x,factor,factorw,factorv,gerr
 128  real(dp) :: ggsum,ggsumsum,ggrestsum
 129  real(dp) :: gijerr,gijsum,gnorm,ln2shx,qphnrm,relchg,tmp,wovert,thmtol
 130  logical :: part1,part2
 131  character(len=500) :: msg
 132  character(len=fnlen) :: thermal_filename
 133  character(len=fnlen) :: bij_filename
 134  character(len=fnlen) :: vij_filename
 135 !arrays
 136  integer :: symrec(3,3,Crystal%nsym),symrel(3,3,Crystal%nsym)
 137  real(dp) :: symrec_cart(3,3,Crystal%nsym)
 138  integer :: igqpt2(3),ii(6),jj(6),qptrlatt(3,3)
 139  integer,allocatable :: indqpt1(:),nchan2(:)
 140  real(dp) :: gprimd(3,3),qphon(3),rprimd(3,3),tens1(3,3),tens2(3,3)
 141  real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom)
 142  real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom)
 143  real(dp) :: phfrq(3*Crystal%natom)
 144  real(dp),allocatable :: bbij(:,:,:),bij(:,:,:),energy(:),energy0(:),entropy(:)
 145  real(dp),allocatable :: entropy0(:),free(:),free0(:),gdos(:,:),gg(:,:),gg_sum(:,:),gg_rest(:,:)
 146  real(dp),allocatable :: ggij(:,:,:,:),gij(:,:,:,:),gw(:,:),qpt2(:,:),spheat(:)
 147  real(dp),allocatable :: spheat0(:),spqpt2(:,:),wme(:),wtq(:),wtq2(:)
 148  real(dp),allocatable :: wtq_folded(:),vij(:,:,:)
 149  real(dp),allocatable :: phon_dos(:)
 150  logical,allocatable :: wgcnv(:),wgijcnv(:)
 151 
 152 ! *********************************************************************
 153 
 154  ! For the time being, this routine can use only 1 MPI process
 155  if (xmpi_comm_rank(comm) /= master) return
 156 
 157  natom = Crystal%natom
 158  symrel = Crystal%symrel
 159  symrec = Crystal%symrec
 160 
 161  call mkrdim(Ifc%acell,Ifc%rprim,rprimd)
 162  call matr3inv(rprimd,gprimd)
 163  do isym = 1, Crystal%nsym
 164    symrec_cart(:,:,isym) = matmul( gprimd, matmul(dble(symrec(:,:,isym)), transpose(rprimd)) )
 165 ! net result is tens1 = rprimd symrec^T gprimd^T   tens1   gprimd symrec rprimd^T
 166  end do
 167 
 168  thermal_filename=trim(outfilename_radix)//"_THERMO"
 169  if (open_file(thermal_filename, msg, newunit=thermal_unit) /= 0) then
 170    MSG_ERROR(msg)
 171  end if
 172 
 173  write(thermal_unit,*) '#'
 174  write(thermal_unit,*) '#  Thermodynamical quantities calculated by ANADDB'
 175  write(thermal_unit,*) '#'
 176 
 177  bij_filename=trim(outfilename_radix)//"_DEBYE_WALLER"
 178  if (open_file(bij_filename, msg, newunit=bij_unit) /= 0) then
 179    MSG_ERROR(msg)
 180  end if
 181 
 182  vij_filename=trim(outfilename_radix)//"_VELOC_SQUARED"
 183  if (open_file(vij_filename, msg, newunit=vij_unit) /= 0) then
 184    MSG_ERROR(msg)
 185  end if
 186 
 187  thmtol = anaddb_dtset%thmtol
 188  nchan=anaddb_dtset%nchan
 189  ntemper=anaddb_dtset%ntemper
 190  nwchan=anaddb_dtset%nwchan
 191  ngrids=anaddb_dtset%ngrids
 192 
 193  ABI_ALLOCATE(bbij,(6,natom,ntemper))
 194  ABI_ALLOCATE(bij,(6,natom,ntemper))
 195  ABI_ALLOCATE(vij,(6,natom,ntemper))
 196  ABI_ALLOCATE(energy,(ntemper))
 197  ABI_ALLOCATE(energy0,(ntemper))
 198  ABI_ALLOCATE(entropy,(ntemper))
 199  ABI_ALLOCATE(entropy0,(ntemper))
 200  ABI_ALLOCATE(free,(ntemper))
 201  ABI_ALLOCATE(free0,(ntemper))
 202 !Doubling the size (nchan) of gg_sum is needed, because the maximum sum frequency is the double
 203 !the for maximum frequency. However, for the write statement, it is better to double
 204 !also the size of the other arrays ...
 205  ABI_ALLOCATE(gdos,(2*nchan,nwchan))
 206  ABI_ALLOCATE(gg,(2*nchan,nwchan))
 207  ABI_ALLOCATE(gg_sum,(2*nchan,nwchan))
 208  ABI_ALLOCATE(gg_rest,(2*nchan,nwchan))
 209  ABI_ALLOCATE(ggij,(6,natom,nchan,nwchan))
 210  ABI_ALLOCATE(gij,(6,natom,nchan,nwchan))
 211  ABI_ALLOCATE(nchan2,(nwchan))
 212  ABI_ALLOCATE(spheat,(ntemper))
 213  ABI_ALLOCATE(spheat0,(ntemper))
 214  ABI_ALLOCATE(wgcnv,(nwchan))
 215  ABI_ALLOCATE(wgijcnv,(nwchan))
 216  ABI_ALLOCATE(wme,(ntemper))
 217  ABI_ALLOCATE(gw,(nchan,nwchan))
 218 
 219 
 220 !initialize ii and jj arrays
 221  ii(1)=1 ; ii(2)=2 ; ii(3)=3
 222  ii(4)=1 ; ii(5)=1 ; ii(6)=2
 223  jj(1)=1 ; jj(2)=2 ; jj(3)=3
 224  jj(4)=2 ; jj(5)=3 ; jj(6)=3
 225 
 226 !Put zeros in the bins of channels of frequencies
 227  gdos(:,:)=0._dp
 228  gg_sum(:,:)=0._dp
 229  gg_rest(:,:)=0._dp
 230  gij(:,:,:,:)=0._dp
 231 
 232 !Neither part1 nor part2 are assumed converged initially
 233 !None of the channel widths are assumed converged initially
 234 
 235  part1=.false.
 236  part2=.false.
 237 
 238  wgcnv(:)=.false.
 239  wgijcnv(:)=.false.
 240 
 241 !Thermodynamic quantities are put to zero.
 242 !(If exactly zero, initial convergence tests will fail.)
 243 
 244  free0(:)=1.d-05
 245  energy0(:)=1.d-05
 246  entropy0(:)=1.d-05
 247  spheat0(:)=1.d-05
 248  bij(:,:,:)=1.d-05
 249  vij(:,:,:)=1.d-05
 250 
 251 !Number of actual channels set
 252  do iwchan=1,nwchan
 253    nchan2(iwchan)=(nchan-1)/iwchan+1
 254  end do
 255 
 256 !For different types of Bravais lattices
 257  facbrv=1
 258  if(anaddb_dtset%brav==2)facbrv=2
 259  if(anaddb_dtset%brav==3)facbrv=4
 260 
 261 !Loops on the q point grids
 262  do igrid=1,ngrids
 263 
 264    igqpt2(:)=max((anaddb_dtset%ng2qpt(:)*igrid)/ngrids, 1)
 265    mqpt2=(igqpt2(1)*igqpt2(2)*igqpt2(3))/facbrv
 266    ABI_ALLOCATE(qpt2,(3,mqpt2))
 267    ABI_ALLOCATE(spqpt2,(3,mqpt2))
 268 
 269    option=1
 270    qptrlatt(:,:)=0
 271    qptrlatt(1,1)=igqpt2(1)
 272    qptrlatt(2,2)=igqpt2(2)
 273    qptrlatt(3,3)=igqpt2(3)
 274    call smpbz(anaddb_dtset%brav,iout,qptrlatt,mqpt2,nspqpt,1,option,anaddb_dtset%q2shft,spqpt2)
 275 
 276    ABI_ALLOCATE(indqpt1,(nspqpt))
 277    ABI_ALLOCATE(wtq,(nspqpt))
 278    ABI_ALLOCATE(wtq_folded,(nspqpt))
 279 
 280 !  Reduce the number of such points by symmetrization
 281    wtq(:)=1.0_dp
 282 
 283    timrev=1
 284    call symkpt(0,Crystal%gmet,indqpt1,ab_out,spqpt2,nspqpt,nqpt2,Crystal%nsym,symrec,timrev,wtq,wtq_folded)
 285    ABI_ALLOCATE(wtq2,(nqpt2))
 286    do iqpt2=1,nqpt2
 287      wtq2(iqpt2)=wtq_folded(indqpt1(iqpt2))
 288      qpt2(:,iqpt2)=spqpt2(:,indqpt1(iqpt2))
 289      !write(std_out,*)' harmonic_thermo : iqpt2, wtq2 :',iqpt2,wtq2(iqpt2)
 290    end do
 291    ABI_DEALLOCATE(wtq_folded)
 292 
 293 !  Temporary counters are put zero.
 294 
 295    gg(:,:)=zero
 296    ggij(:,:,:,:)=zero
 297    gw(:,:)=zero
 298 
 299 !  Sum over the sampled q points
 300    do iqpt2=1,nqpt2
 301 
 302      qphon(:)=qpt2(:,iqpt2)
 303      qphnrm=1.0_dp
 304 
 305      ! Fourier Interpolation
 306      call ifc_fourq(Ifc,Crystal,qphon,phfrq,displ,out_eigvec=eigvec)
 307 
 308      if (present(thmflag)) then
 309        if (thmflag ==2)then
 310          call sortph(eigvec,displ,outfilename_radix,natom,phfrq)
 311        end if
 312      end if
 313 
 314 !    Sum over the phonon modes
 315      do iii=1,3*natom
 316 
 317 !      Slightly negative frequencies are put to zero
 318 !      Imaginary frequencies are also put to zero
 319        if(phfrq(iii)<0._dp) phfrq(iii)=0._dp
 320 
 321 !      Note: frequencies are now in cm-1
 322 
 323 !      Sort frequencies into channels of frequencies for each channel width of frequency
 324        do iwchan=nwchan,1,-1
 325          if (.not.wgcnv(iwchan))then
 326            dosinc=dble(iwchan)
 327            ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
 328 
 329            if(ichan>nchan2(iwchan)) then
 330              write(msg, '(a,es16.6,a,a,a,i7,a,a,a)' )&
 331 &             'There is a phonon frequency,',phfrq(iii),' larger than the maximum one,',ch10,&
 332 &             'as defined by the number of channels of width 1 cm-1, nchan=',nchan,'.',ch10,&
 333 &             'Action: increase nchan (suggestion : double it).'
 334              MSG_ERROR(msg)
 335            end if
 336 
 337            gg(ichan,iwchan)=gg(ichan,iwchan)+wtq2(iqpt2)
 338 
 339            gw(ichan,iwchan)=gw(ichan,iwchan)+wtq2(iqpt2)*phfrq(iii)*Ha_cmm1
 340 
 341 !          to calculate two phonon DOS for qshift = 0.0
 342            do iiii=1,3*natom
 343              if(phfrq(iiii)<0.d0) phfrq(iiii)=0.d0
 344 
 345              ichan=int(abs(phfrq(iii)+phfrq(iiii))*Ha_cmm1/dosinc)+1
 346              gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)+wtq2(iqpt2)
 347 
 348              ichan=int(abs(phfrq(iii)-phfrq(iiii))*Ha_cmm1/dosinc)+1
 349              gg_rest(ichan,iwchan)=gg_rest(ichan,iwchan)+wtq2(iqpt2)
 350            end do ! end loop for iiii
 351 
 352          end if
 353        end do
 354 
 355        do iwchan=nwchan,1,-1
 356          if (.not.wgijcnv(iwchan))then
 357 
 358            dosinc=dble(iwchan)
 359            ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
 360 
 361            if(ichan>nchan2(iwchan)) then
 362              write(msg, '(a,a,a,a,a)' )&
 363 &             'Phonon frequencies larger than the maximum one,',ch10,&
 364 &             'as defined by the number of channels of width 1 cm-1.',ch10,&
 365 &             'Action: increase nchan (suggestion : double it).'
 366              MSG_ERROR(msg)
 367            end if
 368 
 369            do iatom=1,natom
 370              do ij=1,6
 371                ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)&
 372 &               +wtq2(iqpt2)*&
 373 &               (eigvec(1,ii(ij),iatom,iii)*eigvec(1,jj(ij),iatom,iii)+&
 374 &               eigvec(2,ii(ij),iatom,iii)*eigvec(2,jj(ij),iatom,iii) )
 375              end do
 376            end do
 377 
 378          end if
 379        end do
 380 
 381      end do ! End of the sum over the phonon modes
 382    end do ! End of the sum over q-points in the irreducible zone
 383 
 384 !  deallocate sortph array
 385    call end_sortph()
 386 
 387 !  Symmetrize the gij
 388    do ichan=1,nchan
 389      do iwchan=nwchan,1,-1
 390        do iatom=1,natom
 391          do ij=1,6
 392 !          Uses bbij as work space
 393            bbij(ij,iatom,1)=ggij(ij,iatom,ichan,iwchan)
 394            ggij(ij,iatom,ichan,iwchan)=0.0_dp
 395          end do
 396        end do
 397        do iatom=1,natom
 398          do isym=1,Crystal%nsym
 399 !          Find the atom that will be applied on atom iatom
 400            ind=Crystal%indsym(4,isym,iatom)
 401            do ij=1,6
 402              tens1(ii(ij),jj(ij))=bbij(ij,ind,1)
 403            end do
 404 !          complete the 3x3 tensor from the upper triangle
 405            tens1(2,1)=tens1(1,2)
 406            tens1(3,1)=tens1(1,3)
 407            tens1(3,2)=tens1(2,3)
 408 !          Here acomplishes the tensorial operations
 409 !!          make this a BLAS call, or better yet batch the whole thing?
 410 !          2) Apply the symmetry operation on both indices   USING symrec in
 411 !          cartesian coordinates
 412            do iii=1,3
 413              do jjj=1,3
 414                tens2(iii,jjj)=tens1(iii,1)*symrec_cart(1,jjj,isym)&
 415 &               +tens1(iii,2)*symrec_cart(2,jjj,isym)&
 416 &               +tens1(iii,3)*symrec_cart(3,jjj,isym)
 417              end do
 418            end do
 419            do jjj=1,3
 420              do iii=1,3
 421                tens1(iii,jjj)=tens2(1,jjj)*symrec_cart(1,iii,isym)&
 422 &               +tens2(2,jjj)*symrec_cart(2,iii,isym)&
 423 &               +tens2(3,jjj)*symrec_cart(3,iii,isym)
 424              end do
 425            end do
 426 ! net result is tens1 = rprimd symrec^T gprimd^T   tens1   gprimd symrec rprimd^T
 427 
 428 !          This accumulates over atoms, to account for all symmetric ones
 429            do ij=1,6
 430              ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan) + tens1(ii(ij),jj(ij))
 431            end do
 432 
 433          end do
 434 !        Each one will be replicated nsym times in the end:
 435          do ij=1,6
 436            ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)/dble(Crystal%nsym)
 437          end do
 438        end do
 439      end do
 440    end do
 441 
 442    call wrtout(std_out,' harmonic_thermo: g(w) and gij(k|w) calculated given a q sampling grid.','COLL')
 443 
 444 !  Sum up the counts in the channels to check the normalization
 445 !  and test convergence with respect to q sampling
 446    gnorm=(igqpt2(1)*igqpt2(2)*igqpt2(3)*3*natom)/facbrv
 447 
 448    if(.not.part1)then
 449      do iwchan=nwchan,1,-1
 450 
 451        !write(msg,'(a,i0)' )' harmonic_thermo : iwchan=',iwchan
 452        !call wrtout(std_out,msg,'COLL')
 453 
 454        if (wgcnv(iwchan)) cycle
 455        !call wrtout(std_out,' harmonic_thermo : compute g, f, e, s, c','COLL')
 456 
 457 !      Calculate g(w) and F,E,S,C
 458 
 459        ggsum=0.0_dp
 460        ggsumsum=0.0_dp
 461        ggrestsum=0.0_dp
 462        do ichan=1,nchan2(iwchan)
 463          ggsum=ggsum+gg(ichan,iwchan)
 464          ggsumsum=ggsumsum+gg_sum(ichan,iwchan)
 465          ggrestsum=ggrestsum+gg_rest(ichan,iwchan)
 466        end do
 467 
 468        if(ggsum/=gnorm)then
 469          write(msg, '(a,es14.6,i6,a,a,es14.6,a)' )&
 470 &         'Frequencies are missing in g(w) : ggsum,iwchan=',ggsum,iwchan,ch10,&
 471 &         'gnorm=',gnorm,'.'
 472          MSG_BUG(msg)
 473        end if
 474 
 475 !      Check if the density of states changed by more than dostol
 476 
 477        gerr=zero
 478        if (ngrids>1) then
 479          do ichan=1,nchan2(iwchan)
 480            gerr=gerr+abs(gg(ichan,iwchan)/ggsum-gdos(ichan,iwchan))
 481          end do
 482        end if
 483 
 484        if(gerr>anaddb_dtset%dostol.and.ngrids>1) then
 485          wgcnv(iwchan)=.false.
 486        else
 487          wgcnv(iwchan)=.true.
 488        end if
 489 
 490 !      g(w) is updated
 491        do ichan=1,nchan2(iwchan)
 492          gdos(ichan,iwchan)=gg(ichan,iwchan)/ggsum
 493          gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)/ggsumsum
 494          gg_rest(ichan,iwchan)=gg_rest(ichan,iwchan)/ggrestsum
 495        end do
 496        do ichan=1,nchan2(iwchan)
 497          gw(ichan,iwchan)=gw(ichan,iwchan)/ggsum
 498        end do
 499 
 500 !      Write gerr for each q sampling and w width
 501        write(msg,'(a,a,i3,3i6,f10.1,f10.5)') ch10, &
 502 &       'harmonic_thermo: iwchan,igqpt(:),norm,error=',iwchan,igqpt2(1),igqpt2(2),igqpt2(3),ggsum+tol8,gerr+tol10
 503        call wrtout(std_out,msg,'COLL')
 504 
 505 !      If the DOS with a channel width is newly converged,
 506 !      print it out and calculate the thermodynamic functions.
 507        convth=0
 508        if(wgcnv(iwchan)) then
 509          if (ngrids==1) then
 510            if (anaddb_dtset%dossum /= 0 ) then
 511              write(msg,'(a65,i5,a16)') ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS with channel width= ',iwchan,':'
 512            else
 513              write(msg,'(a25,i5,a16)') ' DOS  with channel width=  ',iwchan,':'
 514            end if
 515          else
 516            if (anaddb_dtset%dossum /= 0 ) then
 517              write(msg,'(a65,i5,a16)')&
 518 &             ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS  with channel width= ',iwchan,' newly converged'
 519            else
 520              write(msg,'(a25,i5,a16)') ' DOS  with channel width=  ',iwchan,' newly converged'
 521            end if
 522          end if
 523 
 524          call wrtout(std_out,msg,'COLL')
 525          call wrtout(iout,msg,'COLL')
 526          do ichan=1,nchan2(iwchan)
 527            if (anaddb_dtset%dossum /= 0 ) then
 528              write(msg,'(i8,f11.1,3(f12.5,3x))') ichan,gg(ichan,iwchan)+tol10,&
 529 &             gdos(ichan,iwchan)+tol10, gg_sum(ichan,iwchan)*(3.0*natom*3.0*natom)+tol10, &
 530 &             gg_rest(ichan,iwchan)*(3.0*natom*(3.0*natom-1))+tol10
 531            else
 532              write(msg,'(i8,f11.1,1x,f12.5)') ichan,gg(ichan,iwchan)+tol10,&
 533 &             gdos(ichan,iwchan)+tol10
 534            end if
 535            call wrtout(std_out,msg,'COLL')
 536          end do
 537 
 538          if (ngrids>1) then
 539            write(msg,'(a24,f10.5)')'   with maximal error = ',gerr+tol10
 540            call wrtout(std_out,msg,'COLL')
 541            call wrtout(iout,msg,'COLL')
 542          end if
 543 
 544          nomega = nchan2(iwchan)
 545          dosinc=dble(iwchan)
 546 
 547          ABI_ALLOCATE(phon_dos,(nomega))
 548          phon_dos = gdos(1:nomega,iwchan)
 549 
 550 !Put zeroes for F, E, S, Cv
 551          free(:)=zero
 552          energy(:)=zero
 553          entropy(:)=zero
 554          spheat(:)=zero
 555          wme(:)=zero
 556 
 557          do itemper=1,ntemper
 558 
 559            tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
 560 !          The temperature (tmp) is given in Kelvin
 561            if (tmp < tol6) cycle
 562 
 563            do iomega=1,nomega
 564 
 565 !            wovert= hbar*w / 2kT dimensionless
 566              wovert=dosinc*(dble(iomega)-0.5_dp)/Ha_cmm1/(2._dp*kb_HaK*tmp)
 567              expm2x=exp(-2.0_dp*wovert)
 568              ln2shx=wovert+log(1.0_dp-expm2x)
 569              cothx=(1.0_dp+expm2x)/(1.0_dp-expm2x)
 570              factor=dble(3*natom)*phon_dos(iomega)
 571              factorw=3*natom*gw(iomega,iwchan)
 572 
 573 !            This matches the equations published in Lee & Gonze, PRB 51, 8610 (1995) [[cite:Lee1995]]
 574              free(itemper)=free(itemper) +factor*kb_HaK*tmp*ln2shx
 575              energy(itemper)=energy(itemper) + factor*kb_HaK*tmp*wovert*cothx
 576              entropy(itemper)=entropy(itemper) + factor*kb_HaK*(wovert*cothx - ln2shx)
 577 
 578 !            The contribution is much lower than 1.0d-16 when wovert<100.0_dp
 579              if(wovert<100.0_dp)then
 580                spheat(itemper)=spheat(itemper)+factor*kb_HaK*wovert**2/sinh(wovert)**2
 581              end if
 582              wme(itemper)=wme(itemper)+factorw*kb_HaK*wovert**2/sinh(wovert)**2
 583 
 584            end do ! iomega
 585 
 586            if (abs(spheat(itemper))>tol8) wme(itemper)=wme(itemper)/spheat(itemper)
 587          end do ! itemper
 588          ABI_DEALLOCATE(phon_dos)
 589 
 590 !        Check if the thermodynamic functions change within tolerance,
 591          if (ngrids>1) then
 592            write(msg,'(a,a,a)')&
 593 &           ' harmonic_thermo : check if the thermodynamic functions',ch10,&
 594 &           '    change within tolerance.'
 595            call wrtout(std_out,msg,'COLL')
 596            convth=1
 597 
 598            do itemper=1,ntemper
 599              change=free(itemper)-free0(itemper)
 600              relchg=change/free0(itemper)
 601              if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
 602                write(msg,'(a,es14.4,a,a,es14.4)' )&
 603 &               ' harmonic_thermo : free energy relative changes ',relchg,ch10,&
 604 &               '        are larger than thmtol ',thmtol
 605                call wrtout(std_out,msg,'COLL')
 606                convth=0
 607              end if
 608              change=energy(itemper)-energy0(itemper)
 609              relchg=change/energy0(itemper)
 610              if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
 611                write(msg,'(a,es14.4,a,a,es14.4)' )&
 612 &               ' harmonic_thermo : energy relative changes ',relchg,ch10,&
 613 &               '        are larger than thmtol ',thmtol
 614                call wrtout(std_out,msg,'COLL')
 615                convth=0
 616              end if
 617              change=entropy(itemper)-entropy0(itemper)
 618              relchg=change/entropy0(itemper)
 619              if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
 620                write(msg,'(a,es14.4,a,a,es14.4)' )&
 621 &               ' harmonic_thermo : entropy relative changes ',relchg,ch10,&
 622 &               '        are larger than thmtol ',thmtol
 623                call wrtout(std_out,msg,'COLL')
 624                convth=0
 625              end if
 626              change=spheat(itemper)-spheat0(itemper)
 627              relchg=change/spheat0(itemper)
 628              if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
 629                write(msg,'(a,es14.4,a,a,es14.4)' )&
 630 &               ' harmonic_thermo : specific heat relative changes ',relchg,ch10,&
 631 &               '        are larger than thmtol ',thmtol
 632                call wrtout(std_out,msg,'COLL')
 633                convth=0
 634              end if
 635 
 636              if(convth==0)exit
 637            end do ! End of check if the thermodynamic functions change within tolerance
 638 
 639          else
 640            convth=1
 641          end if
 642 
 643 !        Update F,E,S,C and eventually write them if converged
 644          if(convth==1)then
 645            part1=.true.
 646            write(msg,'(a,a,a)') ch10,&
 647 &           ' # At  T     F(J/mol-c)     E(J/mol-c)     S(J/(mol-c.K)) C(J/(mol-c.K)) Omega_mean(cm-1)'
 648            call wrtout(iout,msg,'COLL')
 649            call wrtout(thermal_unit,msg,'COLL')
 650            msg = ' # (A mol-c is the abbreviation of a mole-cell, that is, the'
 651            call wrtout(iout,msg,'COLL')
 652            call wrtout(thermal_unit,msg,'COLL')
 653            msg = ' #  number of Avogadro times the atoms in a unit cell)'
 654            call wrtout(iout,msg,'COLL')
 655            call wrtout(thermal_unit,msg,'COLL')
 656 
 657            write(msg, '(a,a,a)' )&
 658 &           ' harmonic_thermo : thermodynamic functions have converged',ch10,&
 659 &           '     see main output file ...'
 660            call wrtout(std_out,msg,'COLL')
 661          end if
 662 
 663          do itemper=1,ntemper
 664            free0(itemper)=free(itemper)
 665            energy0(itemper)=energy(itemper)
 666            entropy0(itemper)=entropy(itemper)
 667            spheat0(itemper)=spheat(itemper)
 668 
 669            if(convth==1)then
 670              tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
 671              write(msg,'(es11.3,5es15.7)') tmp+tol8,&
 672 &             Ha_eV*e_Cb*Avogadro*free(itemper),&
 673 &             Ha_eV*e_Cb*Avogadro*energy(itemper),&
 674 &             Ha_eV*e_Cb*Avogadro*entropy(itemper),&
 675 &             Ha_eV*e_Cb*Avogadro*spheat(itemper),&
 676 &             wme(itemper)
 677              call wrtout(iout,msg,'COLL')
 678              call wrtout(thermal_unit,msg,'COLL')
 679            end if
 680          end do
 681        end if
 682 
 683        if(convth==1)exit
 684      end do
 685    end if
 686 
 687    if(.not.part2)then
 688      ! Atomic temperature factor calculation
 689      do iwchan=nwchan,1,-1
 690        if (wgijcnv(iwchan))cycle
 691 
 692        ! Calculate gij(k|w) and Bij(k)
 693        ! Check if the density of states changed by more than dostol
 694        gijsum =zero
 695        wgijcnv(iwchan)=.true.
 696        if (ngrids>1) then
 697          do iatom=1,natom
 698            do ij=1,6
 699              gijerr=zero
 700              do ichan=1,nchan2(iwchan)
 701                gijsum = gijsum + gij(ij,iatom,ichan,iwchan)
 702                gijerr=gijerr&
 703 &               +abs(ggij(ij,iatom,ichan,iwchan)/gnorm&
 704 &               -     gij(ij,iatom,ichan,iwchan))
 705              end do
 706              if(gijerr>anaddb_dtset%dostol) then
 707                wgijcnv(iwchan)=.false.
 708                exit
 709              end if
 710            end do
 711          end do
 712        else
 713          gijerr=0.d0
 714        end if
 715 
 716 !      gij(k|w) is updated
 717 
 718        do ichan=1,nchan2(iwchan)
 719          do iatom=1,natom
 720            do ij=1,6
 721              gij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)/(gnorm/(3*natom))
 722            end do
 723 !if (iwchan==1) write (200+iatom,'(I6,6(E20.10,2x))') ichan, gij(1:6,iatom,ichan,iwchan)
 724          end do
 725        end do
 726 
 727 !      Write gijerr for each q sampling and w width
 728 
 729        write(msg,'(a,a,i3,3i6,f10.5,f10.5)') ch10,&
 730 &       ' iwchan,igqpt(i),gijsum, gij error= ',&
 731 &       iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijsum,gijerr+tol10
 732        call wrtout(std_out,msg,'COLL')
 733 
 734 !      If the generalized DOS with a channel width is newly converged,
 735 !      print it out and calculate Bij(k).
 736        if(wgijcnv(iwchan)) then
 737 
 738          if (ngrids==1) then
 739            write(msg,'(a,i5,a)') ' gij with channel width=  ',iwchan,':'
 740          else
 741            write(msg,'(a,i5,a)') ' gij with channel width=  ',iwchan,' newly converged'
 742          end if
 743          call wrtout(iout,msg,'COLL')
 744 
 745          write(msg,'(a,2i3,3i6,f10.5)')'iatom,iwchan,igqpt2(i),gij error= ',&
 746 &         iatom,iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijerr+tol10
 747          call wrtout(iout,msg,'COLL')
 748 
 749          do itemper=1,ntemper
 750 
 751 !          Put zeroes for Bij(k)
 752            do iatom=1,natom
 753              do ij=1,6
 754                bbij(ij,iatom,itemper)=0._dp
 755                vij(ij,iatom,itemper)=0._dp
 756              end do
 757            end do
 758 
 759            tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
 760 !          tmp in K
 761            if (tmp < tol6) cycle
 762 
 763            dosinc=dble(iwchan)
 764 !
 765            do ichan=1,nchan2(iwchan)
 766 !
 767 !$wovert= \hbar*w / 2kT$, dimensionless
 768              wovert=dosinc*(dble(ichan)-half)/Ha_cmm1/(two*kb_HaK*tmp)
 769              expm2x=exp(-two*wovert)
 770              do iatom=1,natom
 771 !   factor contains 1 / (2 omega)
 772                factor=Ha_cmm1/(two*dosinc*(dble(ichan)-half))    &
 773 &               *(one+expm2x)/(one-expm2x) /amu(Crystal%typat(iatom))/amu_emass
 774 
 775 !   this becomes * 0.5 * omega for the velocities
 776                factorv=(half*dosinc*(dble(ichan)-half)/Ha_cmm1)    &
 777 &               *(one+expm2x)/(one-expm2x) /amu(Crystal%typat(iatom))/amu_emass
 778 
 779                do ij=1,6
 780                  bbij(ij,iatom,itemper)=bbij(ij,iatom,itemper) + factor*gij(ij,iatom,ichan,iwchan)
 781                  vij(ij,iatom,itemper)=vij(ij,iatom,itemper) + factorv*gij(ij,iatom,ichan,iwchan)
 782                end do
 783              end do
 784 
 785            end do
 786 
 787          end do
 788 
 789 !        B matrix is now in atomic unit in the Cartesian coordinates.
 790 !        Check if Bij(k) changed within tolerance.
 791          convth=1
 792          if (ngrids>1) then
 793            do itemper=1,ntemper
 794              do iatom=1,natom
 795                do ij=1,6
 796                  diffbb=bbij(ij,iatom,itemper)-bij(ij,iatom,itemper)
 797                  if (diffbb > 1d-10  .and. diffbb/bij(ij,iatom,itemper) > thmtol) then
 798                    write(msg,'(a)' )' harmonic_thermo : Bij changes are larger than thmtol '
 799                    call wrtout(std_out,msg,'COLL')
 800                    convth=0
 801                  end if
 802                  if(convth==0)exit
 803                end do
 804                if(convth==0)exit
 805              end do
 806              if(convth==0)exit
 807            end do
 808          end if
 809 
 810          bij=bbij ! save for next iteration
 811 
 812          !Update Bij(k) and write them. B matrix printed in angstrom^2
 813          !TODO : get rid of this version in the log and output file. Prefer
 814          !external files
 815          if (convth==1) then
 816            write(msg, '(a,a,a)' )&
 817 &           ' B matrix elements as a function of T',ch10,&
 818 &           '    Angstrom^2, cartesian coordinates'
 819            call wrtout(std_out,msg,'COLL')
 820            call wrtout(iout,msg,'COLL')
 821 
 822            do itemper=1,ntemper
 823 !            tmp in K
 824              tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
 825              do iatom=1,natom
 826                write(iout,'(2i3,es11.3,6es12.4)')&
 827 &               iwchan,iatom,tmp+tol10,&
 828 &               Bohr_Ang**2*bij(1,iatom,itemper)+tol10,&
 829 &               Bohr_Ang**2*bij(2,iatom,itemper)+tol10,&
 830 &               Bohr_Ang**2*bij(3,iatom,itemper)+tol10,&
 831 &               Bohr_Ang**2*bij(4,iatom,itemper)+tol10,&
 832 &               Bohr_Ang**2*bij(5,iatom,itemper)+tol10,&
 833 &               Bohr_Ang**2*bij(6,iatom,itemper)+tol10
 834              end do ! end loop over natom
 835            end do ! end loop over ntemper
 836 
 837 !        Mean square velocity matrix printed in angstrom^2/picosec^2
 838            write(msg, '(a,a,a)' )&
 839 &           ' <vel^2> matrix elements as a function of T',ch10,&
 840 &           '    Angstrom^2/(picosec)^2, cartesian coordinates'
 841            call wrtout(std_out,msg,'COLL')
 842            call wrtout(iout,msg,'COLL')
 843 
 844            do itemper=1,ntemper
 845 !            tmp in K
 846              tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
 847              do iatom=1,natom
 848                vij(:,iatom,itemper)=Bohr_Ang**2*vij(:,iatom,itemper)/(Time_Sec*1.0e12)**2
 849 !              The following check zeros out <v^2> if it is very small, in order to
 850 !              avoid numerical noise being interpreted by the automatic tests as
 851 !              something real. Note also that we compare it in
 852 !              absolute value, that's because if any of the phonon frequencies are
 853 !              computed as negative, <v^2> can take a negative value.
 854                do icomp=1, 6
 855                  if (abs(vij(icomp,iatom,itemper)) < 1.0e-12) vij(icomp,iatom,itemper)=zero
 856                end do
 857                write(iout,'(2i3,es11.3,6es12.4)')&
 858 &               iwchan,iatom,tmp+tol10,&
 859 &               vij(1,iatom,itemper),&
 860 &               vij(2,iatom,itemper),&
 861 &               vij(3,iatom,itemper),&
 862 &               vij(4,iatom,itemper),&
 863 &               vij(5,iatom,itemper),&
 864 &               vij(6,iatom,itemper)
 865              end do ! end loop over natom
 866            end do ! end loop over ntemper
 867          end if ! end check on convergence
 868 
 869 
 870          ! keep this one !!!!!!!!!!!!!!!!!!
 871          if (convth==1) then
 872            write(msg, '(a,a,a)' )&
 873 &           '# B matrix elements as a function of T, for each atom, and smallest omega channel width',ch10,&
 874 &           '#    Angstrom^2, cartesian coordinates'
 875            call wrtout(bij_unit,msg,'COLL')
 876            do iatom=1,natom
 877              write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
 878              call wrtout(bij_unit,msg,'COLL')
 879              do itemper=1,ntemper
 880 !              tmp in K
 881                tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
 882                write(msg,'(es11.3,6es12.4)')&
 883 &               tmp,&
 884 &               Bohr_Ang**2*bij(1,iatom,itemper),&
 885 &               Bohr_Ang**2*bij(2,iatom,itemper),&
 886 &               Bohr_Ang**2*bij(3,iatom,itemper),&
 887 &               Bohr_Ang**2*bij(4,iatom,itemper),&
 888 &               Bohr_Ang**2*bij(5,iatom,itemper),&
 889 &               Bohr_Ang**2*bij(6,iatom,itemper)
 890                call wrtout(bij_unit,msg,'COLL')
 891              end do ! end loop over ntemper
 892            end do ! end loop over natom
 893 
 894 !        Mean square velocity matrix printed in angstrom^2/picosec^2
 895            write(msg, '(a,a,a)' )&
 896 &           '# <vel^2> matrix elements as a function of T, for each atom, and smallest channel width',ch10,&
 897 &           '#    Angstrom^2/(picosec)^2, cartesian coordinates'
 898            call wrtout(vij_unit,msg,'COLL')
 899 
 900            do iatom=1,natom
 901              write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
 902              call wrtout(vij_unit,msg,'COLL')
 903              do itemper=1,ntemper
 904 !            tmp in K
 905                tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
 906                vij(:,iatom,itemper)=Bohr_Ang**2*vij(:,iatom,itemper)/(Time_Sec*1.0e12)**2
 907 
 908 !            The following check zeros out <v^2> if it is very small, in order to
 909 !            avoid numerical noise being interpreted by the automatic tests as
 910 !            something real. Note also that we compare it in
 911 !            absolute value, that's because if any of the phonon frequencies are
 912 !            computed as negative, <v^2> can take a negative value.
 913                do icomp=1, 6
 914                  if (abs(vij(icomp,iatom,itemper)) < 1.0e-12) vij(icomp,iatom,itemper)=zero
 915                end do
 916                write(vij_unit,'(es11.3,6es12.4)')&
 917 &               tmp,&
 918 &               vij(1,iatom,itemper),&
 919 &               vij(2,iatom,itemper),&
 920 &               vij(3,iatom,itemper),&
 921 &               vij(4,iatom,itemper),&
 922 &               vij(5,iatom,itemper),&
 923 &               vij(6,iatom,itemper)
 924              end do ! end loop over ntemper
 925            end do ! end loop over natom
 926          end if ! end check on convergence
 927 
 928          if(convth==1)part2=.true.
 929 
 930        end if ! End of test on wgijcnv
 931      end do ! End of loop over iwchan
 932    end if ! End of part2
 933 
 934    if(part1.and.part2)exit
 935 
 936    ABI_DEALLOCATE(indqpt1)
 937    ABI_DEALLOCATE(qpt2)
 938    ABI_DEALLOCATE(spqpt2)
 939    ABI_DEALLOCATE(wtq)
 940    ABI_DEALLOCATE(wtq2)
 941 
 942  end do ! End of the Loop on the q point grids
 943 
 944  ABI_DEALLOCATE(bbij)
 945  ABI_DEALLOCATE(bij)
 946  ABI_DEALLOCATE(energy)
 947  ABI_DEALLOCATE(energy0)
 948  ABI_DEALLOCATE(entropy)
 949  ABI_DEALLOCATE(entropy0)
 950  ABI_DEALLOCATE(free)
 951  ABI_DEALLOCATE(free0)
 952  ABI_DEALLOCATE(gdos)
 953  ABI_DEALLOCATE(gg)
 954  ABI_DEALLOCATE(gg_sum)
 955  ABI_DEALLOCATE(gg_rest)
 956  ABI_DEALLOCATE(ggij)
 957  ABI_DEALLOCATE(gij)
 958  ABI_DEALLOCATE(nchan2)
 959  ABI_DEALLOCATE(spheat)
 960  ABI_DEALLOCATE(spheat0)
 961  ABI_DEALLOCATE(vij)
 962  ABI_DEALLOCATE(wgcnv)
 963  ABI_DEALLOCATE(wgijcnv)
 964  if(allocated(indqpt1)) then
 965    ABI_DEALLOCATE(indqpt1)
 966  end if
 967  if(allocated(qpt2)) then
 968    ABI_DEALLOCATE(qpt2)
 969  end if
 970  if(allocated(spqpt2)) then
 971    ABI_DEALLOCATE(spqpt2)
 972  end if
 973  if(allocated(wtq)) then
 974    ABI_DEALLOCATE(wtq)
 975  end if
 976  if(allocated(wtq2)) then
 977    ABI_DEALLOCATE(wtq2)
 978  end if
 979  ABI_DEALLOCATE(gw)
 980  ABI_DEALLOCATE(wme)
 981 
 982  if(.not.part1)then
 983    write(msg, '(a,a,a,a,a,a,a,a,a)' )&
 984 &   'No thermodynamical function is printed out :',ch10,&
 985 &   'the tolerance level that was asked ',ch10,&
 986 &   'has not been match with the grids specified.',ch10,&
 987 &   'Action: in the input file, increase the resolution',ch10,&
 988 &   'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
 989    MSG_ERROR(msg)
 990  end if
 991 
 992  if(.not.part2)then
 993    write(msg,'(a,a,a,a,a,a,a,a,a)')&
 994 &   'No atomic factor tensor is printed out :',ch10,&
 995 &   'the tolerance level that was asked ',ch10,&
 996 &   'has not been match with the grids specified.',ch10,&
 997 &   'Action: in the input file, increase the resolution',ch10,&
 998 &   'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
 999    MSG_WARNING(msg)
1000  end if
1001 
1002  close (thermal_unit)
1003  close (bij_unit)
1004  close (vij_unit)
1005 
1006 end subroutine harmonic_thermo