TABLE OF CONTENTS
ABINIT/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
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