Branch data 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 : : !!***
|
