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