TABLE OF CONTENTS


ABINIT/harmonic_thermo [ Functions ]

[ Top ] [ Functions ]

NAME

 harmonic_thermo

FUNCTION

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

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (CL,XG)
 This file is distributed under the terms of the
 GNU General Public Licence, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

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