Line data Source code
1 : !!****m* ABINIT/m_harmonic_thermo
2 : !! NAME
3 : !! m_harmonic_thermo
4 : !!
5 : !! FUNCTION
6 : !! This routine to calculate phonon density of states,
7 : !! thermodynamical properties, Debye-Waller factor, and atomic mean square velocity
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (CL, XG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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_matrix, 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
44 : !!***
45 :
46 : public :: harmonic_thermo
47 : !!***
48 :
49 : contains
50 : !!***
51 :
52 : !!****f* m_harmonic_thermo/harmonic_thermo
53 : !!
54 : !! NAME
55 : !! harmonic_thermo
56 : !!
57 : !! FUNCTION
58 : !! This routine to calculate phonon density of states,
59 : !! thermodynamical properties, Debye-Waller factor, and atomic mean square velocity
60 : !!
61 : !! INPUTS
62 : !! Crystal<crystal_t>=data type gathering info on the crystalline structure.
63 : !! Ifc<ifc_type>=Object containing the interatomic force constants.
64 : !! %atmfrc(2,3,natom,3,natom,nrpt) = Interatomic Forces in real space
65 : !! %dyewq0(3,3,natom)=Ewald part of the dynamical matrix, at q=0.
66 : !! %rpt(3,nrpt)=canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!)
67 : !! %nrpt=number of R points in the Big Box
68 : !! %trans(3,natom)=atomic translations : xred = rcan + trans
69 : !! %wghatm(natom,natom,nrpt)=weights associated to a pair of atoms and to a R vector
70 : !! amu(ntypat)=mass of the atoms (atomic mass unit)
71 : !! anaddb_dtset= (derived datatype) contains all the input variables
72 : !! iout =unit number for output
73 : !! natom=number of atoms in the unit cell
74 : !! outfilename_radix=radix of anaddb output file name: append _THERMO for thermodynamic quantities
75 : !! comm=MPI communicator
76 : !!
77 : !! OUTPUT
78 : !!
79 : !! NOTES
80 : !! dosinc=increment between the channels for the phonon DOS in cm-1
81 : !!
82 : !! SOURCE
83 :
84 4 : 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 8 : integer :: symrec(3,3,Crystal%nsym),symrel(3,3,Crystal%nsym)
117 8 : real(dp) :: symrec_cart(3,3,Crystal%nsym)
118 : integer :: igqpt2(3),ii(6),jj(6),qptrlatt(3,3)
119 4 : integer,allocatable :: indqpt1(:),nchan2(:),bz2ibz_smap(:,:)
120 : real(dp) :: gprimd(3,3),qphon(3),rprimd(3,3),tens1(3,3),tens2(3,3)
121 8 : real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom)
122 8 : real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom)
123 4 : real(dp) :: phfrq(3*Crystal%natom)
124 4 : real(dp),allocatable :: bbij(:,:,:),bij(:,:,:),energy(:),energy0(:),entropy(:)
125 4 : real(dp),allocatable :: entropy0(:),free(:),free0(:),gdos(:,:),gg(:,:),gg_sum(:,:),gg_rest(:,:)
126 4 : real(dp),allocatable :: ggij(:,:,:,:),gij(:,:,:,:),gw(:,:),qpt2(:,:),spheat(:)
127 4 : real(dp),allocatable :: spheat0(:),spqpt2(:,:),wme(:),wtq(:),wtq2(:)
128 4 : real(dp),allocatable :: wtq_folded(:),vij(:,:,:)
129 4 : real(dp),allocatable :: phon_dos(:)
130 4 : logical,allocatable :: wgcnv(:),wgijcnv(:)
131 :
132 : ! *********************************************************************
133 :
134 : ! For the time being, this routine can use only 1 MPI process
135 4 : if (xmpi_comm_rank(comm) /= master) return
136 :
137 4 : natom = Crystal%natom
138 784 : symrel = Crystal%symrel
139 784 : symrec = Crystal%symrec
140 :
141 4 : call mkrdim(Ifc%acell,Ifc%rprim,rprimd)
142 4 : call matr3inv(rprimd,gprimd)
143 64 : do isym = 1, Crystal%nsym
144 6184 : 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 4 : thermal_filename=trim(outfilename_radix)//"_THERMO"
149 4 : if (open_file(thermal_filename, msg, newunit=thermal_unit) /= 0) then
150 0 : ABI_ERROR(msg)
151 : end if
152 :
153 4 : write(thermal_unit,*) '#'
154 4 : write(thermal_unit,*) '# Thermodynamical quantities calculated by ANADDB'
155 4 : write(thermal_unit,*) '#'
156 :
157 4 : bij_filename=trim(outfilename_radix)//"_DEBYE_WALLER"
158 4 : if (open_file(bij_filename, msg, newunit=bij_unit) /= 0) then
159 0 : ABI_ERROR(msg)
160 : end if
161 :
162 4 : vij_filename=trim(outfilename_radix)//"_VELOC_SQUARED"
163 4 : if (open_file(vij_filename, msg, newunit=vij_unit) /= 0) then
164 0 : ABI_ERROR(msg)
165 : end if
166 :
167 4 : thmtol = anaddb_dtset%thmtol
168 4 : nchan=anaddb_dtset%nchan
169 4 : ntemper=anaddb_dtset%ntemper
170 4 : nwchan=anaddb_dtset%nwchan
171 4 : ngrids=anaddb_dtset%ngrids
172 :
173 16 : ABI_MALLOC(bbij,(6,natom,ntemper))
174 12 : ABI_MALLOC(bij,(6,natom,ntemper))
175 12 : ABI_MALLOC(vij,(6,natom,ntemper))
176 12 : ABI_MALLOC(energy,(ntemper))
177 8 : ABI_MALLOC(energy0,(ntemper))
178 8 : ABI_MALLOC(entropy,(ntemper))
179 8 : ABI_MALLOC(entropy0,(ntemper))
180 8 : ABI_MALLOC(free,(ntemper))
181 8 : 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 16 : ABI_MALLOC(gdos,(2*nchan,nwchan))
186 12 : ABI_MALLOC(gg,(2*nchan,nwchan))
187 12 : ABI_MALLOC(gg_sum,(2*nchan,nwchan))
188 12 : ABI_MALLOC(gg_rest,(2*nchan,nwchan))
189 20 : ABI_MALLOC(ggij,(6,natom,nchan,nwchan))
190 16 : ABI_MALLOC(gij,(6,natom,nchan,nwchan))
191 12 : ABI_MALLOC(nchan2,(nwchan))
192 8 : ABI_MALLOC(spheat,(ntemper))
193 8 : ABI_MALLOC(spheat0,(ntemper))
194 8 : ABI_MALLOC(wgcnv,(nwchan))
195 8 : ABI_MALLOC(wgijcnv,(nwchan))
196 8 : ABI_MALLOC(wme,(ntemper))
197 16 : ABI_MALLOC(gw,(nchan,nwchan))
198 :
199 :
200 : !initialize ii and jj arrays
201 4 : ii(1)=1 ; ii(2)=2 ; ii(3)=3
202 4 : ii(4)=1 ; ii(5)=1 ; ii(6)=2
203 4 : jj(1)=1 ; jj(2)=2 ; jj(3)=3
204 4 : jj(4)=2 ; jj(5)=3 ; jj(6)=3
205 :
206 : !Put zeros in the bins of channels of frequencies
207 53529 : gdos(:,:)=0._dp
208 53529 : gg_sum(:,:)=0._dp
209 53529 : gg_rest(:,:)=0._dp
210 1013779 : 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 29 : part1=.false.
216 29 : part2=.false.
217 :
218 29 : wgcnv(:)=.false.
219 29 : wgijcnv(:)=.false.
220 :
221 : !Thermodynamic quantities are put to zero.
222 : !(If exactly zero, initial convergence tests will fail.)
223 :
224 44 : free0(:)=1.d-05
225 44 : energy0(:)=1.d-05
226 44 : entropy0(:)=1.d-05
227 44 : spheat0(:)=1.d-05
228 1584 : bij(:,:,:)=1.d-05
229 1584 : vij(:,:,:)=1.d-05
230 :
231 : !Number of actual channels set
232 29 : do iwchan=1,nwchan
233 29 : nchan2(iwchan)=(nchan-1)/iwchan+1
234 : end do
235 :
236 : !For different types of Bravais lattices
237 4 : facbrv=1
238 4 : if(anaddb_dtset%brav==2)facbrv=2
239 4 : if(anaddb_dtset%brav==3)facbrv=4
240 :
241 : !Loops on the q point grids
242 16 : do igrid=1,ngrids
243 :
244 56 : igqpt2(:)=max((anaddb_dtset%ng2qpt(:)*igrid)/ngrids, 1)
245 14 : mqpt2=(igqpt2(1)*igqpt2(2)*igqpt2(3))/facbrv
246 42 : ABI_MALLOC(qpt2,(3,mqpt2))
247 28 : ABI_MALLOC(spqpt2,(3,mqpt2))
248 :
249 14 : option=1
250 14 : qptrlatt(:,:)=0
251 14 : qptrlatt(1,1)=igqpt2(1)
252 14 : qptrlatt(2,2)=igqpt2(2)
253 14 : qptrlatt(3,3)=igqpt2(3)
254 14 : call smpbz(anaddb_dtset%brav,iout,qptrlatt,mqpt2,nspqpt,1,option,anaddb_dtset%q2shft,spqpt2)
255 :
256 42 : ABI_MALLOC(indqpt1,(nspqpt))
257 42 : ABI_MALLOC(wtq,(nspqpt))
258 28 : ABI_MALLOC(wtq_folded,(nspqpt))
259 42 : ABI_MALLOC(bz2ibz_smap,(6, nspqpt))
260 :
261 : ! Reduce the number of such points by symmetrization
262 31918 : wtq(:)=1.0_dp
263 :
264 14 : timrev=1
265 : call symkpt(0,Crystal%gmet,indqpt1,ab_out,spqpt2,nspqpt,nqpt2,Crystal%nsym,symrec,timrev,wtq,wtq_folded, &
266 14 : bz2ibz_smap, xmpi_comm_self)
267 :
268 14 : ABI_FREE(bz2ibz_smap)
269 :
270 42 : ABI_MALLOC(wtq2,(nqpt2))
271 1760 : do iqpt2=1,nqpt2
272 1746 : wtq2(iqpt2)=wtq_folded(indqpt1(iqpt2))
273 6998 : qpt2(:,iqpt2)=spqpt2(:,indqpt1(iqpt2))
274 : !write(std_out,*)' harmonic_thermo : iqpt2, wtq2 :',iqpt2,wtq2(iqpt2)
275 : end do
276 14 : ABI_FREE(wtq_folded)
277 :
278 : ! Temporary counters are put zero.
279 :
280 192609 : gg(:,:)=zero
281 2668859 : ggij(:,:,:,:)=zero
282 96359 : gw(:,:)=zero
283 :
284 : ! Sum over the sampled q points
285 1760 : do iqpt2=1,nqpt2
286 :
287 6984 : qphon(:)=qpt2(:,iqpt2)
288 1746 : qphnrm=1.0_dp
289 :
290 : ! Fourier Interpolation
291 1746 : call ifc%fourq(Crystal,qphon,phfrq,displ,out_eigvec=eigvec)
292 :
293 1746 : if (present(thmflag)) then
294 0 : if (thmflag ==2)then
295 0 : call sortph(eigvec,displ,outfilename_radix,natom,phfrq)
296 : end if
297 : end if
298 :
299 : ! Sum over the phonon modes
300 33782 : do iii=1,3*natom
301 :
302 : ! Slightly negative frequencies are put to zero
303 : ! Imaginary frequencies are also put to zero
304 32022 : 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 198432 : do iwchan=nwchan,1,-1
310 198432 : if (.not.wgcnv(iwchan))then
311 164874 : dosinc=dble(iwchan)
312 164874 : ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
313 :
314 164874 : if(ichan>nchan2(iwchan)) then
315 : write(msg, '(a,es16.6,a,a,a,i7,a,a,a)' )&
316 0 : & 'There is a phonon frequency,',phfrq(iii),' larger than the maximum one,',ch10,&
317 0 : & 'as defined by the number of channels of width 1 cm-1, nchan=',nchan,'.',ch10,&
318 0 : & 'Action: increase nchan (suggestion : double it).'
319 0 : ABI_ERROR(msg)
320 : end if
321 :
322 164874 : gg(ichan,iwchan)=gg(ichan,iwchan)+wtq2(iqpt2)
323 :
324 164874 : gw(ichan,iwchan)=gw(ichan,iwchan)+wtq2(iqpt2)*phfrq(iii)*Ha_cmm1
325 :
326 : ! to calculate two phonon DOS for qshift = 0.0
327 4062828 : do iiii=1,3*natom
328 3897954 : if(phfrq(iiii)<0.d0) phfrq(iiii)=0.d0
329 :
330 3897954 : ichan=int(abs(phfrq(iii)+phfrq(iiii))*Ha_cmm1/dosinc)+1
331 3897954 : gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)+wtq2(iqpt2)
332 :
333 3897954 : ichan=int(abs(phfrq(iii)-phfrq(iiii))*Ha_cmm1/dosinc)+1
334 4062828 : 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 200178 : do iwchan=nwchan,1,-1
341 198432 : if (.not.wgijcnv(iwchan))then
342 :
343 105660 : dosinc=dble(iwchan)
344 105660 : ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
345 :
346 105660 : if(ichan>nchan2(iwchan)) then
347 : write(msg, '(a,a,a,a,a)' )&
348 0 : & 'Phonon frequencies larger than the maximum one,',ch10,&
349 0 : & 'as defined by the number of channels of width 1 cm-1.',ch10,&
350 0 : & 'Action: increase nchan (suggestion : double it).'
351 0 : ABI_ERROR(msg)
352 : end if
353 :
354 1054080 : do iatom=1,natom
355 6744600 : 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 6638940 : & 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 14 : call end_sortph()
371 :
372 : ! Symmetrize the gij
373 15264 : do ichan=1,nchan
374 111514 : do iwchan=nwchan,1,-1
375 463750 : do iatom=1,natom
376 2668750 : do ij=1,6
377 : ! Uses bbij as work space
378 2205000 : bbij(ij,iatom,1)=ggij(ij,iatom,ichan,iwchan)
379 2572500 : ggij(ij,iatom,ichan,iwchan)=0.0_dp
380 : end do
381 : end do
382 479000 : do iatom=1,natom
383 5137500 : do isym=1,Crystal%nsym
384 : ! Find the atom that will be applied on atom iatom
385 4770000 : ind=Crystal%indsym(4,isym,iatom)
386 33390000 : do ij=1,6
387 33390000 : tens1(ii(ij),jj(ij))=bbij(ij,ind,1)
388 : end do
389 : ! complete the 3x3 tensor from the upper triangle
390 4770000 : tens1(2,1)=tens1(1,2)
391 4770000 : tens1(3,1)=tens1(1,3)
392 4770000 : 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 19080000 : do iii=1,3
398 62010000 : 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 57240000 : & +tens1(iii,3)*symrec_cart(3,jjj,isym)
402 : end do
403 : end do
404 19080000 : do jjj=1,3
405 62010000 : 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 57240000 : & +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 33757500 : do ij=1,6
415 33390000 : 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 2668750 : do ij=1,6
421 2572500 : 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 14 : 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 14 : gnorm=(igqpt2(1)*igqpt2(2)*igqpt2(3)*3*natom)/facbrv
432 :
433 14 : if(.not.part1)then
434 87 : do iwchan=nwchan,1,-1
435 :
436 : !write(msg,'(a,i0)' )' harmonic_thermo : iwchan=',iwchan
437 : !call wrtout(std_out,msg,'COLL')
438 :
439 77 : 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 76 : ggsum=0.0_dp
445 76 : ggsumsum=0.0_dp
446 76 : ggrestsum=0.0_dp
447 27881 : do ichan=1,nchan2(iwchan)
448 27805 : ggsum=ggsum+gg(ichan,iwchan)
449 27805 : ggsumsum=ggsumsum+gg_sum(ichan,iwchan)
450 27881 : ggrestsum=ggrestsum+gg_rest(ichan,iwchan)
451 : end do
452 :
453 76 : if(ggsum/=gnorm)then
454 : write(msg, '(a,es14.6,i6,a,a,es14.6,a)' )&
455 0 : & 'Frequencies are missing in g(w) : ggsum,iwchan=',ggsum,iwchan,ch10,&
456 0 : & 'gnorm=',gnorm,'.'
457 0 : ABI_BUG(msg)
458 : end if
459 :
460 : ! Check if the density of states changed by more than dostol
461 :
462 76 : gerr=zero
463 76 : if (ngrids>1) then
464 27630 : do ichan=1,nchan2(iwchan)
465 27630 : gerr=gerr+abs(gg(ichan,iwchan)/ggsum-gdos(ichan,iwchan))
466 : end do
467 : end if
468 :
469 76 : if(gerr>anaddb_dtset%dostol.and.ngrids>1) then
470 69 : wgcnv(iwchan)=.false.
471 : else
472 7 : wgcnv(iwchan)=.true.
473 : end if
474 :
475 : ! g(w) is updated
476 27881 : do ichan=1,nchan2(iwchan)
477 27805 : gdos(ichan,iwchan)=gg(ichan,iwchan)/ggsum
478 27805 : gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)/ggsumsum
479 27881 : gg_rest(ichan,iwchan)=gg_rest(ichan,iwchan)/ggrestsum
480 : end do
481 27881 : do ichan=1,nchan2(iwchan)
482 27881 : gw(ichan,iwchan)=gw(ichan,iwchan)/ggsum
483 : end do
484 :
485 : ! Write gerr for each q sampling and w width
486 76 : write(msg,'(a,a,i3,3i6,f10.1,f10.5)') ch10, &
487 152 : & 'harmonic_thermo: iwchan,igqpt(:),norm,error=',iwchan,igqpt2(1),igqpt2(2),igqpt2(3),ggsum+tol8,gerr+tol10
488 76 : 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 76 : convth=0
493 76 : if(wgcnv(iwchan)) then
494 7 : if (ngrids==1) then
495 1 : if (anaddb_dtset%dossum /= 0 ) then
496 0 : write(msg,'(a65,i5,a16)') ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS with channel width= ',iwchan,':'
497 : else
498 1 : write(msg,'(a25,i5,a16)') ' DOS with channel width= ',iwchan,':'
499 : end if
500 : else
501 6 : if (anaddb_dtset%dossum /= 0 ) then
502 : write(msg,'(a65,i5,a16)')&
503 2 : & ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS with channel width= ',iwchan,' newly converged'
504 : else
505 4 : write(msg,'(a25,i5,a16)') ' DOS with channel width= ',iwchan,' newly converged'
506 : end if
507 : end if
508 :
509 7 : call wrtout(std_out,msg,'COLL')
510 7 : call wrtout(iout,msg,'COLL')
511 1552 : do ichan=1,nchan2(iwchan)
512 1545 : if (anaddb_dtset%dossum /= 0 ) then
513 563 : write(msg,'(i8,f11.1,3(f12.5,3x))') ichan,gg(ichan,iwchan)+tol10,&
514 563 : & gdos(ichan,iwchan)+tol10, gg_sum(ichan,iwchan)*(3.0*natom*3.0*natom)+tol10, &
515 1126 : & gg_rest(ichan,iwchan)*(3.0*natom*(3.0*natom-1))+tol10
516 : else
517 982 : write(msg,'(i8,f11.1,1x,f12.5)') ichan,gg(ichan,iwchan)+tol10,&
518 1964 : & gdos(ichan,iwchan)+tol10
519 : end if
520 1552 : call wrtout(std_out,msg,'COLL')
521 : end do
522 :
523 7 : if (ngrids>1) then
524 6 : write(msg,'(a24,f10.5)')' with maximal error = ',gerr+tol10
525 6 : call wrtout(std_out,msg,'COLL')
526 6 : call wrtout(iout,msg,'COLL')
527 : end if
528 :
529 7 : nomega = nchan2(iwchan)
530 7 : dosinc=dble(iwchan)
531 :
532 21 : ABI_MALLOC(phon_dos,(nomega))
533 1559 : phon_dos = gdos(1:nomega,iwchan)
534 :
535 : !Put zeroes for F, E, S, Cv
536 77 : free(:)=zero
537 77 : energy(:)=zero
538 77 : entropy(:)=zero
539 77 : spheat(:)=zero
540 77 : wme(:)=zero
541 :
542 77 : do itemper=1,ntemper
543 :
544 70 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
545 : ! The temperature (tmp) is given in Kelvin
546 70 : if (tmp < tol6) cycle
547 :
548 15520 : do iomega=1,nomega
549 :
550 : ! wovert= hbar*w / 2kT dimensionless
551 15450 : wovert=dosinc*(dble(iomega)-0.5_dp)/Ha_cmm1/(2._dp*kb_HaK*tmp)
552 15450 : expm2x=exp(-2.0_dp*wovert)
553 15450 : ln2shx=wovert+log(1.0_dp-expm2x)
554 15450 : cothx=(1.0_dp+expm2x)/(1.0_dp-expm2x)
555 15450 : factor=dble(3*natom)*phon_dos(iomega)
556 15450 : factorw=3*natom*gw(iomega,iwchan)
557 :
558 : ! This matches the equations published in Lee & Gonze, PRB 51, 8610 (1995) [[cite:Lee1995]]
559 15450 : free(itemper)=free(itemper) +factor*kb_HaK*tmp*ln2shx
560 15450 : energy(itemper)=energy(itemper) + factor*kb_HaK*tmp*wovert*cothx
561 15450 : 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 15450 : if(wovert<100.0_dp)then
565 15450 : spheat(itemper)=spheat(itemper)+factor*kb_HaK*wovert**2/sinh(wovert)**2
566 : end if
567 15520 : wme(itemper)=wme(itemper)+factorw*kb_HaK*wovert**2/sinh(wovert)**2
568 :
569 : end do ! iomega
570 :
571 77 : if (abs(spheat(itemper))>tol8) wme(itemper)=wme(itemper)/spheat(itemper)
572 : end do ! itemper
573 7 : ABI_FREE(phon_dos)
574 :
575 : ! Check if the thermodynamic functions change within tolerance,
576 7 : if (ngrids>1) then
577 : write(msg,'(a,a,a)')&
578 6 : & ' harmonic_thermo : check if the thermodynamic functions',ch10,&
579 12 : & ' change within tolerance.'
580 6 : call wrtout(std_out,msg,'COLL')
581 6 : convth=1
582 :
583 36 : do itemper=1,ntemper
584 33 : change=free(itemper)-free0(itemper)
585 33 : relchg=change/free0(itemper)
586 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
587 : write(msg,'(a,es14.4,a,a,es14.4)' )&
588 3 : & ' harmonic_thermo : free energy relative changes ',relchg,ch10,&
589 6 : & ' are larger than thmtol ',thmtol
590 3 : call wrtout(std_out,msg,'COLL')
591 3 : convth=0
592 : end if
593 33 : change=energy(itemper)-energy0(itemper)
594 33 : relchg=change/energy0(itemper)
595 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
596 : write(msg,'(a,es14.4,a,a,es14.4)' )&
597 3 : & ' harmonic_thermo : energy relative changes ',relchg,ch10,&
598 6 : & ' are larger than thmtol ',thmtol
599 3 : call wrtout(std_out,msg,'COLL')
600 3 : convth=0
601 : end if
602 33 : change=entropy(itemper)-entropy0(itemper)
603 33 : relchg=change/entropy0(itemper)
604 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
605 : write(msg,'(a,es14.4,a,a,es14.4)' )&
606 1 : & ' harmonic_thermo : entropy relative changes ',relchg,ch10,&
607 2 : & ' are larger than thmtol ',thmtol
608 1 : call wrtout(std_out,msg,'COLL')
609 1 : convth=0
610 : end if
611 33 : change=spheat(itemper)-spheat0(itemper)
612 33 : relchg=change/spheat0(itemper)
613 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
614 : write(msg,'(a,es14.4,a,a,es14.4)' )&
615 1 : & ' harmonic_thermo : specific heat relative changes ',relchg,ch10,&
616 2 : & ' are larger than thmtol ',thmtol
617 1 : call wrtout(std_out,msg,'COLL')
618 : convth=0
619 : end if
620 :
621 35 : 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 6 : if(convth==1)then
630 4 : part1=.true.
631 4 : write(msg,'(a,a,a)') ch10,&
632 8 : & ' # 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 4 : call wrtout(iout,msg,'COLL')
634 4 : call wrtout(thermal_unit,msg,'COLL')
635 4 : msg = ' # (A mol-c is the abbreviation of a mole-cell, that is, the'
636 4 : call wrtout(iout,msg,'COLL')
637 4 : call wrtout(thermal_unit,msg,'COLL')
638 4 : msg = ' # number of Avogadro times the atoms in a unit cell)'
639 4 : call wrtout(iout,msg,'COLL')
640 4 : call wrtout(thermal_unit,msg,'COLL')
641 :
642 : write(msg, '(a,a,a)' )&
643 4 : & ' harmonic_thermo : thermodynamic functions have converged',ch10,&
644 8 : & ' see main output file ...'
645 4 : call wrtout(std_out,msg,'COLL')
646 : end if
647 :
648 77 : do itemper=1,ntemper
649 70 : free0(itemper)=free(itemper)
650 70 : energy0(itemper)=energy(itemper)
651 70 : entropy0(itemper)=entropy(itemper)
652 70 : spheat0(itemper)=spheat(itemper)
653 :
654 77 : if(convth==1)then
655 40 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
656 40 : write(msg,'(es11.3,5es15.7)') tmp+tol8,&
657 40 : & Ha_eV*e_Cb*Avogadro*free(itemper),&
658 40 : & Ha_eV*e_Cb*Avogadro*energy(itemper),&
659 40 : & Ha_eV*e_Cb*Avogadro*entropy(itemper),&
660 40 : & Ha_eV*e_Cb*Avogadro*spheat(itemper),&
661 80 : & wme(itemper)
662 40 : call wrtout(iout,msg,'COLL')
663 40 : call wrtout(thermal_unit,msg,'COLL')
664 : end if
665 : end do
666 : end if
667 :
668 17 : if(convth==1)exit
669 : end do
670 : end if
671 :
672 14 : if(.not.part2)then
673 : ! Atomic temperature factor calculation
674 65 : do iwchan=nwchan,1,-1
675 55 : 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 25 : gijsum =zero
680 25 : wgijcnv(iwchan)=.true.
681 25 : if (ngrids>1) then
682 95 : do iatom=1,natom
683 545 : do ij=1,6
684 450 : gijerr=zero
685 217020 : do ichan=1,nchan2(iwchan)
686 216570 : gijsum = gijsum + gij(ij,iatom,ichan,iwchan)
687 : gijerr=gijerr&
688 : & +abs(ggij(ij,iatom,ichan,iwchan)/gnorm&
689 217020 : & - gij(ij,iatom,ichan,iwchan))
690 : end do
691 525 : if(gijerr>anaddb_dtset%dostol) then
692 0 : wgijcnv(iwchan)=.false.
693 0 : 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 10935 : do ichan=1,nchan2(iwchan)
704 72725 : do iatom=1,natom
705 443440 : do ij=1,6
706 432530 : 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 25 : write(msg,'(a,a,i3,3i6,f10.5,f10.5)') ch10,&
715 25 : & ' iwchan,igqpt(i),gijsum, gij error= ',&
716 50 : & iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijsum,gijerr+tol10
717 25 : 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 35 : if(wgijcnv(iwchan)) then
722 :
723 25 : if (ngrids==1) then
724 5 : write(msg,'(a,i5,a)') ' gij with channel width= ',iwchan,':'
725 : else
726 20 : write(msg,'(a,i5,a)') ' gij with channel width= ',iwchan,' newly converged'
727 : end if
728 25 : call wrtout(iout,msg,'COLL')
729 :
730 25 : write(msg,'(a,2i3,3i6,f10.5)')'iatom,iwchan,igqpt2(i),gij error= ',&
731 50 : & iatom,iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijerr+tol10
732 25 : call wrtout(iout,msg,'COLL')
733 :
734 275 : do itemper=1,ntemper
735 :
736 : ! Put zeroes for Bij(k)
737 1450 : do iatom=1,natom
738 8650 : do ij=1,6
739 7200 : bbij(ij,iatom,itemper)=0._dp
740 8400 : vij(ij,iatom,itemper)=0._dp
741 : end do
742 : end do
743 :
744 250 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
745 : ! tmp in K
746 250 : if (tmp < tol6) cycle
747 :
748 250 : dosinc=dble(iwchan)
749 : !
750 109375 : do ichan=1,nchan2(iwchan)
751 : !
752 : !$wovert= \hbar*w / 2kT$, dimensionless
753 109100 : wovert=dosinc*(dble(ichan)-half)/Ha_cmm1/(two*kb_HaK*tmp)
754 109100 : expm2x=exp(-two*wovert)
755 727250 : do iatom=1,natom
756 : ! factor contains 1 / (2 omega)
757 : factor=Ha_cmm1/(two*dosinc*(dble(ichan)-half)) &
758 617900 : & *(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 617900 : & *(one+expm2x)/(one-expm2x) /amu(Crystal%typat(iatom))/amu_emass
763 :
764 4434400 : do ij=1,6
765 3707400 : bbij(ij,iatom,itemper)=bbij(ij,iatom,itemper) + factor*gij(ij,iatom,ichan,iwchan)
766 4325300 : 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 25 : convth=1
777 25 : if (ngrids>1) then
778 112 : do itemper=1,ntemper
779 309 : do iatom=1,natom
780 1467 : do ij=1,6
781 1261 : 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 1261 : if (diffbb > 1d-10) then
788 409 : if (bij(ij,iatom,itemper) /= zero) then
789 409 : if (diffbb/bij(ij,iatom,itemper) > thmtol) then
790 11 : write(msg,'(a)' )' harmonic_thermo : Bij changes are larger than thmtol '
791 11 : call wrtout(std_out,msg,'COLL')
792 : convth=0
793 : end if
794 : end if
795 : end if
796 206 : if(convth==0)exit
797 : end do
798 309 : if(convth==0)exit
799 : end do
800 112 : if(convth==0)exit
801 : end do
802 : end if
803 :
804 8700 : 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 25 : if (convth==1) then
810 : write(msg, '(a,a,a)' )&
811 14 : & ' B matrix elements as a function of T',ch10,&
812 28 : & ' Angstrom^2, cartesian coordinates'
813 14 : call wrtout(std_out,msg,'COLL')
814 14 : call wrtout(iout,msg,'COLL')
815 :
816 154 : do itemper=1,ntemper
817 : ! tmp in K
818 140 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
819 784 : do iatom=1,natom
820 : write(iout,'(2i3,es11.3,6es12.4)')&
821 630 : & iwchan,iatom,tmp+tol10,&
822 630 : & Bohr_Ang**2*bij(1,iatom,itemper)+tol10,&
823 630 : & Bohr_Ang**2*bij(2,iatom,itemper)+tol10,&
824 630 : & Bohr_Ang**2*bij(3,iatom,itemper)+tol10,&
825 630 : & Bohr_Ang**2*bij(4,iatom,itemper)+tol10,&
826 630 : & Bohr_Ang**2*bij(5,iatom,itemper)+tol10,&
827 1400 : & 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 14 : & ' <vel^2> matrix elements as a function of T',ch10,&
834 28 : & ' Angstrom^2/(picosec)^2, cartesian coordinates'
835 14 : call wrtout(std_out,msg,'COLL')
836 14 : call wrtout(iout,msg,'COLL')
837 :
838 154 : do itemper=1,ntemper
839 : ! tmp in K
840 140 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
841 784 : do iatom=1,natom
842 4410 : 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 4410 : do icomp=1, 6
849 4410 : 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 630 : & iwchan,iatom,tmp+tol10,&
853 630 : & vij(1,iatom,itemper),&
854 630 : & vij(2,iatom,itemper),&
855 630 : & vij(3,iatom,itemper),&
856 630 : & vij(4,iatom,itemper),&
857 630 : & vij(5,iatom,itemper),&
858 1400 : & 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 14 : & '# B matrix elements as a function of T, for each atom, and smallest omega channel width',ch10,&
868 28 : & '# Angstrom^2, cartesian coordinates'
869 14 : call wrtout(bij_unit,msg,'COLL')
870 77 : do iatom=1,natom
871 63 : write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
872 63 : call wrtout(bij_unit,msg,'COLL')
873 707 : do itemper=1,ntemper
874 : ! tmp in K
875 630 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
876 : write(msg,'(es11.3,6es12.4)')&
877 630 : & tmp,&
878 630 : & Bohr_Ang**2*bij(1,iatom,itemper),&
879 630 : & Bohr_Ang**2*bij(2,iatom,itemper),&
880 630 : & Bohr_Ang**2*bij(3,iatom,itemper),&
881 630 : & Bohr_Ang**2*bij(4,iatom,itemper),&
882 630 : & Bohr_Ang**2*bij(5,iatom,itemper),&
883 1260 : & Bohr_Ang**2*bij(6,iatom,itemper)
884 693 : 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 14 : & '# <vel^2> matrix elements as a function of T, for each atom, and smallest channel width',ch10,&
891 28 : & '# Angstrom^2/(picosec)^2, cartesian coordinates'
892 14 : call wrtout(vij_unit,msg,'COLL')
893 :
894 77 : do iatom=1,natom
895 63 : write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
896 63 : call wrtout(vij_unit,msg,'COLL')
897 707 : do itemper=1,ntemper
898 : ! tmp in K
899 630 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
900 4410 : 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 4410 : do icomp=1, 6
908 4410 : 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 630 : & tmp,&
912 630 : & vij(1,iatom,itemper),&
913 630 : & vij(2,iatom,itemper),&
914 630 : & vij(3,iatom,itemper),&
915 630 : & vij(4,iatom,itemper),&
916 630 : & vij(5,iatom,itemper),&
917 1323 : & 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 14 : if(part1.and.part2)exit
929 :
930 12 : ABI_FREE(indqpt1)
931 12 : ABI_FREE(qpt2)
932 12 : ABI_FREE(spqpt2)
933 12 : ABI_FREE(wtq)
934 16 : ABI_FREE(wtq2)
935 :
936 : end do ! End of the Loop on the q point grids
937 :
938 4 : ABI_FREE(bbij)
939 4 : ABI_FREE(bij)
940 4 : ABI_FREE(energy)
941 4 : ABI_FREE(energy0)
942 4 : ABI_FREE(entropy)
943 4 : ABI_FREE(entropy0)
944 4 : ABI_FREE(free)
945 4 : ABI_FREE(free0)
946 4 : ABI_FREE(gdos)
947 4 : ABI_FREE(gg)
948 4 : ABI_FREE(gg_sum)
949 4 : ABI_FREE(gg_rest)
950 4 : ABI_FREE(ggij)
951 4 : ABI_FREE(gij)
952 4 : ABI_FREE(nchan2)
953 4 : ABI_FREE(spheat)
954 4 : ABI_FREE(spheat0)
955 4 : ABI_FREE(vij)
956 4 : ABI_FREE(wgcnv)
957 4 : ABI_FREE(wgijcnv)
958 4 : if(allocated(indqpt1)) then
959 2 : ABI_FREE(indqpt1)
960 : end if
961 4 : if(allocated(qpt2)) then
962 2 : ABI_FREE(qpt2)
963 : end if
964 4 : if(allocated(spqpt2)) then
965 2 : ABI_FREE(spqpt2)
966 : end if
967 4 : if(allocated(wtq)) then
968 2 : ABI_FREE(wtq)
969 : end if
970 4 : if(allocated(wtq2)) then
971 2 : ABI_FREE(wtq2)
972 : end if
973 4 : ABI_FREE(gw)
974 4 : ABI_FREE(wme)
975 :
976 4 : if(.not.part1)then
977 : write(msg, '(a,a,a,a,a,a,a,a,a)' )&
978 0 : & 'No thermodynamical function is printed out :',ch10,&
979 0 : & 'the tolerance level that was asked ',ch10,&
980 0 : & 'has not been match with the grids specified.',ch10,&
981 0 : & 'Action: in the input file, increase the resolution',ch10,&
982 0 : & 'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
983 0 : ABI_ERROR(msg)
984 : end if
985 :
986 4 : if(.not.part2)then
987 : write(msg,'(a,a,a,a,a,a,a,a,a)')&
988 2 : & 'No atomic factor tensor is printed out :',ch10,&
989 2 : & 'the tolerance level that was asked ',ch10,&
990 2 : & 'has not been match with the grids specified.',ch10,&
991 2 : & 'Action: in the input file, increase the resolution',ch10,&
992 4 : & 'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
993 2 : ABI_WARNING(msg)
994 : end if
995 :
996 4 : close (thermal_unit)
997 4 : close (bij_unit)
998 4 : close (vij_unit)
999 :
1000 : end subroutine harmonic_thermo
1001 : !!***
1002 :
1003 : end module m_harmonic_thermo
1004 : !!***
|
