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, qmesh_from_brav_ngqpt, kptrlatt_ncells
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,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,nqshft_use
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),q2shft_use(3,1)
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 [ + + + + ]: 52 : rprimd = Ifc%rprimd
142 [ + + + + ]: 52 : gprimd = Ifc%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 : : !Loops on the q point grids
237 [ + + ]: 16 : do igrid=1,ngrids
238 : :
239 [ + + ]: 56 : igqpt2(:)=max((anaddb_dtset%ng2qpt(:)*igrid)/ngrids, 1)
240 : : ! Translate the (brav, ngqpt, q2shft) specification into a brav-independent q-mesh lattice.
241 : 14 : option=1
242 : 14 : call qmesh_from_brav_ngqpt(anaddb_dtset%brav, igqpt2, 1, anaddb_dtset%q2shft, qptrlatt, nqshft_use, q2shft_use)
243 : 14 : mqpt2 = kptrlatt_ncells(qptrlatt, nqshft_use)
244 [ - + + - : 42 : ABI_MALLOC(qpt2,(3,mqpt2))
+ - + - -
+ - + ]
245 [ - + + - : 28 : ABI_MALLOC(spqpt2,(3,mqpt2))
- + - + ]
246 : 14 : call smpbz(iout,qptrlatt,mqpt2,nspqpt,nqshft_use,option,q2shft_use,spqpt2)
247 : :
248 [ - + + - : 42 : ABI_MALLOC(indqpt1,(nspqpt))
+ - + - -
+ - + ]
249 [ - + + - : 42 : ABI_MALLOC(wtq,(nspqpt))
+ - + - -
+ - + ]
250 [ - + + - : 28 : ABI_MALLOC(wtq_folded,(nspqpt))
- + - + ]
251 [ - + + - : 42 : ABI_MALLOC(bz2ibz_smap,(6, nspqpt))
+ - + - -
+ - + ]
252 : :
253 : : ! Reduce the number of such points by symmetrization
254 [ + + ]: 31918 : wtq(:)=1.0_dp
255 : :
256 : 14 : timrev=1
257 : : call symkpt(0,Crystal%gmet,indqpt1,ab_out,spqpt2,nspqpt,nqpt2,Crystal%nsym,symrec,timrev,wtq,wtq_folded, &
258 : 14 : bz2ibz_smap, xmpi_comm_self)
259 : :
260 [ - + ]: 14 : ABI_FREE(bz2ibz_smap)
261 : :
262 [ - + + - : 42 : ABI_MALLOC(wtq2,(nqpt2))
+ - + - -
+ - + ]
263 [ + + ]: 1760 : do iqpt2=1,nqpt2
264 : 1746 : wtq2(iqpt2)=wtq_folded(indqpt1(iqpt2))
265 [ + + ]: 6998 : qpt2(:,iqpt2)=spqpt2(:,indqpt1(iqpt2))
266 : : !write(std_out,*)' harmonic_thermo : iqpt2, wtq2 :',iqpt2,wtq2(iqpt2)
267 : : end do
268 [ - + ]: 14 : ABI_FREE(wtq_folded)
269 : :
270 : : ! Temporary counters are put zero.
271 : :
272 [ + + + + ]: 192609 : gg(:,:)=zero
273 [ + + + + : 2668859 : ggij(:,:,:,:)=zero
+ + + + ]
274 [ + + + + ]: 96359 : gw(:,:)=zero
275 : :
276 : : ! Sum over the sampled q points
277 [ + + ]: 1760 : do iqpt2=1,nqpt2
278 : :
279 [ + + ]: 6984 : qphon(:)=qpt2(:,iqpt2)
280 : 1746 : qphnrm=1.0_dp
281 : :
282 : : ! Fourier Interpolation
283 : 1746 : call ifc%fourq(Crystal,qphon,phfrq,displ,out_eigvec=eigvec)
284 : :
285 [ - + ]: 1746 : if (present(thmflag)) then
286 [ # # ]: 0 : if (thmflag ==2)then
287 : 0 : call sortph(eigvec,displ,outfilename_radix,natom,phfrq)
288 : : end if
289 : : end if
290 : :
291 : : ! Sum over the phonon modes
292 [ + + ]: 33782 : do iii=1,3*natom
293 : :
294 : : ! Slightly negative frequencies are put to zero
295 : : ! Imaginary frequencies are also put to zero
296 [ + + ]: 32022 : if(phfrq(iii)<0._dp) phfrq(iii)=0._dp
297 : :
298 : : ! Note: frequencies are now in cm-1
299 : :
300 : : ! Sort frequencies into channels of frequencies for each channel width of frequency
301 [ + + ]: 198432 : do iwchan=nwchan,1,-1
302 [ + + ]: 198432 : if (.not.wgcnv(iwchan))then
303 : 164874 : dosinc=dble(iwchan)
304 : 164874 : ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
305 : :
306 [ - + ]: 164874 : if(ichan>nchan2(iwchan)) then
307 : : write(msg, '(a,es16.6,a,a,a,i7,a,a,a)' )&
308 : 0 : & 'There is a phonon frequency,',phfrq(iii),' larger than the maximum one,',ch10,&
309 : 0 : & 'as defined by the number of channels of width 1 cm-1, nchan=',nchan,'.',ch10,&
310 : 0 : & 'Action: increase nchan (suggestion : double it).'
311 : 0 : ABI_ERROR(msg)
312 : : end if
313 : :
314 : 164874 : gg(ichan,iwchan)=gg(ichan,iwchan)+wtq2(iqpt2)
315 : :
316 : 164874 : gw(ichan,iwchan)=gw(ichan,iwchan)+wtq2(iqpt2)*phfrq(iii)*Ha_cmm1
317 : :
318 : : ! to calculate two phonon DOS for qshift = 0.0
319 [ + + ]: 4062828 : do iiii=1,3*natom
320 [ - + ]: 3897954 : if(phfrq(iiii)<0.d0) phfrq(iiii)=0.d0
321 : :
322 : 3897954 : ichan=int(abs(phfrq(iii)+phfrq(iiii))*Ha_cmm1/dosinc)+1
323 : 3897954 : gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)+wtq2(iqpt2)
324 : :
325 : 3897954 : ichan=int(abs(phfrq(iii)-phfrq(iiii))*Ha_cmm1/dosinc)+1
326 : 4062828 : gg_rest(ichan,iwchan)=gg_rest(ichan,iwchan)+wtq2(iqpt2)
327 : : end do ! end loop for iiii
328 : :
329 : : end if
330 : : end do
331 : :
332 [ + + ]: 200178 : do iwchan=nwchan,1,-1
333 [ + + ]: 198432 : if (.not.wgijcnv(iwchan))then
334 : :
335 : 105660 : dosinc=dble(iwchan)
336 : 105660 : ichan=int(phfrq(iii)*Ha_cmm1/dosinc)+1
337 : :
338 [ - + ]: 105660 : if(ichan>nchan2(iwchan)) then
339 : : write(msg, '(a,a,a,a,a)' )&
340 : 0 : & 'Phonon frequencies larger than the maximum one,',ch10,&
341 : 0 : & 'as defined by the number of channels of width 1 cm-1.',ch10,&
342 : 0 : & 'Action: increase nchan (suggestion : double it).'
343 : 0 : ABI_ERROR(msg)
344 : : end if
345 : :
346 [ + + ]: 1054080 : do iatom=1,natom
347 [ + + ]: 6744600 : do ij=1,6
348 : : ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)&
349 : : & +wtq2(iqpt2)*&
350 : : & (eigvec(1,ii(ij),iatom,iii)*eigvec(1,jj(ij),iatom,iii)+&
351 : 6638940 : & eigvec(2,ii(ij),iatom,iii)*eigvec(2,jj(ij),iatom,iii) )
352 : : end do
353 : : end do
354 : :
355 : : end if
356 : : end do
357 : :
358 : : end do ! End of the sum over the phonon modes
359 : : end do ! End of the sum over q-points in the irreducible zone
360 : :
361 : : ! deallocate sortph array
362 : 14 : call end_sortph()
363 : :
364 : : ! Symmetrize the gij
365 [ + + ]: 15264 : do ichan=1,nchan
366 [ + + ]: 111514 : do iwchan=nwchan,1,-1
367 [ + + ]: 463750 : do iatom=1,natom
368 [ + + ]: 2668750 : do ij=1,6
369 : : ! Uses bbij as work space
370 : 2205000 : bbij(ij,iatom,1)=ggij(ij,iatom,ichan,iwchan)
371 : 2572500 : ggij(ij,iatom,ichan,iwchan)=0.0_dp
372 : : end do
373 : : end do
374 [ + + ]: 479000 : do iatom=1,natom
375 [ + + ]: 5137500 : do isym=1,Crystal%nsym
376 : : ! Find the atom that will be applied on atom iatom
377 : 4770000 : ind=Crystal%indsym(4,isym,iatom)
378 [ + + ]: 33390000 : do ij=1,6
379 : 33390000 : tens1(ii(ij),jj(ij))=bbij(ij,ind,1)
380 : : end do
381 : : ! complete the 3x3 tensor from the upper triangle
382 : 4770000 : tens1(2,1)=tens1(1,2)
383 : 4770000 : tens1(3,1)=tens1(1,3)
384 : 4770000 : tens1(3,2)=tens1(2,3)
385 : : ! Here acomplishes the tensorial operations
386 : : !! make this a BLAS call, or better yet batch the whole thing?
387 : : ! 2) Apply the symmetry operation on both indices USING symrec in
388 : : ! cartesian coordinates
389 [ + + ]: 19080000 : do iii=1,3
390 [ + + ]: 62010000 : do jjj=1,3
391 : : tens2(iii,jjj)=tens1(iii,1)*symrec_cart(1,jjj,isym)&
392 : : & +tens1(iii,2)*symrec_cart(2,jjj,isym)&
393 : 57240000 : & +tens1(iii,3)*symrec_cart(3,jjj,isym)
394 : : end do
395 : : end do
396 [ + + ]: 19080000 : do jjj=1,3
397 [ + + ]: 62010000 : do iii=1,3
398 : : tens1(iii,jjj)=tens2(1,jjj)*symrec_cart(1,iii,isym)&
399 : : & +tens2(2,jjj)*symrec_cart(2,iii,isym)&
400 : 57240000 : & +tens2(3,jjj)*symrec_cart(3,iii,isym)
401 : : end do
402 : : end do
403 : : ! net result is tens1 = rprimd symrec^T gprimd^T tens1 gprimd symrec rprimd^T
404 : :
405 : : ! This accumulates over atoms, to account for all symmetric ones
406 [ + + ]: 33757500 : do ij=1,6
407 : 33390000 : ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan) + tens1(ii(ij),jj(ij))
408 : : end do
409 : :
410 : : end do
411 : : ! Each one will be replicated nsym times in the end:
412 [ + + ]: 2668750 : do ij=1,6
413 : 2572500 : ggij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)/dble(Crystal%nsym)
414 : : end do
415 : : end do
416 : : end do
417 : : end do
418 : :
419 : 14 : call wrtout(std_out,' harmonic_thermo: g(w) and gij(k|w) calculated given a q sampling grid.','COLL')
420 : :
421 : : ! Sum up the counts in the channels to check the normalization
422 : : ! and test convergence with respect to q sampling
423 : 14 : gnorm=nspqpt*3*natom
424 : :
425 [ + - ]: 14 : if(.not.part1)then
426 [ + + ]: 87 : do iwchan=nwchan,1,-1
427 : :
428 : : !write(msg,'(a,i0)' )' harmonic_thermo : iwchan=',iwchan
429 : : !call wrtout(std_out,msg,'COLL')
430 : :
431 [ + + ]: 77 : if (wgcnv(iwchan)) cycle
432 : : !call wrtout(std_out,' harmonic_thermo : compute g, f, e, s, c','COLL')
433 : :
434 : : ! Calculate g(w) and F,E,S,C
435 : :
436 : 76 : ggsum=0.0_dp
437 : 76 : ggsumsum=0.0_dp
438 : 76 : ggrestsum=0.0_dp
439 [ + + ]: 27881 : do ichan=1,nchan2(iwchan)
440 : 27805 : ggsum=ggsum+gg(ichan,iwchan)
441 : 27805 : ggsumsum=ggsumsum+gg_sum(ichan,iwchan)
442 : 27881 : ggrestsum=ggrestsum+gg_rest(ichan,iwchan)
443 : : end do
444 : :
445 [ - + ]: 76 : if(ggsum/=gnorm)then
446 : : write(msg, '(a,es14.6,i6,a,a,es14.6,a)' )&
447 : 0 : & 'Frequencies are missing in g(w) : ggsum,iwchan=',ggsum,iwchan,ch10,&
448 : 0 : & 'gnorm=',gnorm,'.'
449 : 0 : ABI_BUG(msg)
450 : : end if
451 : :
452 : : ! Check if the density of states changed by more than dostol
453 : :
454 : 76 : gerr=zero
455 [ + + ]: 76 : if (ngrids>1) then
456 [ + + ]: 27630 : do ichan=1,nchan2(iwchan)
457 : 27630 : gerr=gerr+abs(gg(ichan,iwchan)/ggsum-gdos(ichan,iwchan))
458 : : end do
459 : : end if
460 : :
461 [ + + + - ]: 76 : if(gerr>anaddb_dtset%dostol.and.ngrids>1) then
462 : 69 : wgcnv(iwchan)=.false.
463 : : else
464 : 7 : wgcnv(iwchan)=.true.
465 : : end if
466 : :
467 : : ! g(w) is updated
468 [ + + ]: 27881 : do ichan=1,nchan2(iwchan)
469 : 27805 : gdos(ichan,iwchan)=gg(ichan,iwchan)/ggsum
470 : 27805 : gg_sum(ichan,iwchan)=gg_sum(ichan,iwchan)/ggsumsum
471 : 27881 : gg_rest(ichan,iwchan)=gg_rest(ichan,iwchan)/ggrestsum
472 : : end do
473 [ + + ]: 27881 : do ichan=1,nchan2(iwchan)
474 : 27881 : gw(ichan,iwchan)=gw(ichan,iwchan)/ggsum
475 : : end do
476 : :
477 : : ! Write gerr for each q sampling and w width
478 : 76 : write(msg,'(a,a,i3,3i6,f10.1,f10.5)') ch10, &
479 : 152 : & 'harmonic_thermo: iwchan,igqpt(:),norm,error=',iwchan,igqpt2(1),igqpt2(2),igqpt2(3),ggsum+tol8,gerr+tol10
480 : 76 : call wrtout(std_out,msg,'COLL')
481 : :
482 : : ! If the DOS with a channel width is newly converged,
483 : : ! print it out and calculate the thermodynamic functions.
484 : 76 : convth=0
485 [ + + ]: 76 : if(wgcnv(iwchan)) then
486 [ + + ]: 7 : if (ngrids==1) then
487 [ - + ]: 1 : if (anaddb_dtset%dossum /= 0 ) then
488 : 0 : write(msg,'(a65,i5,a16)') ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS with channel width= ',iwchan,':'
489 : : else
490 : 1 : write(msg,'(a25,i5,a16)') ' DOS with channel width= ',iwchan,':'
491 : : end if
492 : : else
493 [ + + ]: 6 : if (anaddb_dtset%dossum /= 0 ) then
494 : : write(msg,'(a65,i5,a16)')&
495 : 2 : & ' DOS, SUM AND DIFFERENCE OF TWO PHONON DOS with channel width= ',iwchan,' newly converged'
496 : : else
497 : 4 : write(msg,'(a25,i5,a16)') ' DOS with channel width= ',iwchan,' newly converged'
498 : : end if
499 : : end if
500 : :
501 : 7 : call wrtout(std_out,msg,'COLL')
502 : 7 : call wrtout(iout,msg,'COLL')
503 [ + + ]: 1552 : do ichan=1,nchan2(iwchan)
504 [ + + ]: 1545 : if (anaddb_dtset%dossum /= 0 ) then
505 : 563 : write(msg,'(i8,f11.1,3(f12.5,3x))') ichan,gg(ichan,iwchan)+tol10,&
506 : 563 : & gdos(ichan,iwchan)+tol10, gg_sum(ichan,iwchan)*(3.0*natom*3.0*natom)+tol10, &
507 : 1126 : & gg_rest(ichan,iwchan)*(3.0*natom*(3.0*natom-1))+tol10
508 : : else
509 : 982 : write(msg,'(i8,f11.1,1x,f12.5)') ichan,gg(ichan,iwchan)+tol10,&
510 : 1964 : & gdos(ichan,iwchan)+tol10
511 : : end if
512 : 1552 : call wrtout(std_out,msg,'COLL')
513 : : end do
514 : :
515 [ + + ]: 7 : if (ngrids>1) then
516 : 6 : write(msg,'(a24,f10.5)')' with maximal error = ',gerr+tol10
517 : 6 : call wrtout(std_out,msg,'COLL')
518 : 6 : call wrtout(iout,msg,'COLL')
519 : : end if
520 : :
521 : 7 : nomega = nchan2(iwchan)
522 : 7 : dosinc=dble(iwchan)
523 : :
524 [ - + + - : 21 : ABI_MALLOC(phon_dos,(nomega))
+ - + - -
+ - + ]
525 [ - + + + ]: 1559 : phon_dos = gdos(1:nomega,iwchan)
526 : :
527 : : !Put zeroes for F, E, S, Cv
528 [ + + ]: 77 : free(:)=zero
529 [ + + ]: 77 : energy(:)=zero
530 [ + + ]: 77 : entropy(:)=zero
531 [ + + ]: 77 : spheat(:)=zero
532 [ + + ]: 77 : wme(:)=zero
533 : :
534 [ + + ]: 77 : do itemper=1,ntemper
535 : :
536 : 70 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
537 : : ! The temperature (tmp) is given in Kelvin
538 [ + - ]: 70 : if (tmp < tol6) cycle
539 : :
540 [ + + ]: 15520 : do iomega=1,nomega
541 : :
542 : : ! wovert= hbar*w / 2kT dimensionless
543 : 15450 : wovert=dosinc*(dble(iomega)-0.5_dp)/Ha_cmm1/(2._dp*kb_HaK*tmp)
544 : 15450 : expm2x=exp(-2.0_dp*wovert)
545 : 15450 : ln2shx=wovert+log(1.0_dp-expm2x)
546 : 15450 : cothx=(1.0_dp+expm2x)/(1.0_dp-expm2x)
547 : 15450 : factor=dble(3*natom)*phon_dos(iomega)
548 : 15450 : factorw=3*natom*gw(iomega,iwchan)
549 : :
550 : : ! This matches the equations published in Lee & Gonze, PRB 51, 8610 (1995) [[cite:Lee1995]]
551 : 15450 : free(itemper)=free(itemper) +factor*kb_HaK*tmp*ln2shx
552 : 15450 : energy(itemper)=energy(itemper) + factor*kb_HaK*tmp*wovert*cothx
553 : 15450 : entropy(itemper)=entropy(itemper) + factor*kb_HaK*(wovert*cothx - ln2shx)
554 : :
555 : : ! The contribution is much lower than 1.0d-16 when wovert<100.0_dp
556 [ + - ]: 15450 : if(wovert<100.0_dp)then
557 : 15450 : spheat(itemper)=spheat(itemper)+factor*kb_HaK*wovert**2/sinh(wovert)**2
558 : : end if
559 : 15520 : wme(itemper)=wme(itemper)+factorw*kb_HaK*wovert**2/sinh(wovert)**2
560 : :
561 : : end do ! iomega
562 : :
563 [ + - ]: 77 : if (abs(spheat(itemper))>tol8) wme(itemper)=wme(itemper)/spheat(itemper)
564 : : end do ! itemper
565 [ - + ]: 7 : ABI_FREE(phon_dos)
566 : :
567 : : ! Check if the thermodynamic functions change within tolerance,
568 [ + + ]: 7 : if (ngrids>1) then
569 : : write(msg,'(a,a,a)')&
570 : 6 : & ' harmonic_thermo : check if the thermodynamic functions',ch10,&
571 : 12 : & ' change within tolerance.'
572 : 6 : call wrtout(std_out,msg,'COLL')
573 : 6 : convth=1
574 : :
575 [ + + ]: 36 : do itemper=1,ntemper
576 : 33 : change=free(itemper)-free0(itemper)
577 : 33 : relchg=change/free0(itemper)
578 [ + - + + ]: 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
579 : : write(msg,'(a,es14.4,a,a,es14.4)' )&
580 : 3 : & ' harmonic_thermo : free energy relative changes ',relchg,ch10,&
581 : 6 : & ' are larger than thmtol ',thmtol
582 : 3 : call wrtout(std_out,msg,'COLL')
583 : 3 : convth=0
584 : : end if
585 : 33 : change=energy(itemper)-energy0(itemper)
586 : 33 : relchg=change/energy0(itemper)
587 [ + + + + ]: 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
588 : : write(msg,'(a,es14.4,a,a,es14.4)' )&
589 : 3 : & ' harmonic_thermo : energy relative changes ',relchg,ch10,&
590 : 6 : & ' are larger than thmtol ',thmtol
591 : 3 : call wrtout(std_out,msg,'COLL')
592 : 3 : convth=0
593 : : end if
594 : 33 : change=entropy(itemper)-entropy0(itemper)
595 : 33 : relchg=change/entropy0(itemper)
596 [ + + + - ]: 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
597 : : write(msg,'(a,es14.4,a,a,es14.4)' )&
598 : 1 : & ' harmonic_thermo : entropy relative changes ',relchg,ch10,&
599 : 2 : & ' are larger than thmtol ',thmtol
600 : 1 : call wrtout(std_out,msg,'COLL')
601 : 1 : convth=0
602 : : end if
603 : 33 : change=spheat(itemper)-spheat0(itemper)
604 : 33 : relchg=change/spheat0(itemper)
605 [ + + + + ]: 33 : if(change>1d-14*dble(mqpt2) .and. relchg>thmtol)then
606 : : write(msg,'(a,es14.4,a,a,es14.4)' )&
607 : 1 : & ' harmonic_thermo : specific heat relative changes ',relchg,ch10,&
608 : 2 : & ' are larger than thmtol ',thmtol
609 : 1 : call wrtout(std_out,msg,'COLL')
610 : : convth=0
611 : : end if
612 : :
613 [ + + ]: 35 : if(convth==0)exit
614 : : end do ! End of check if the thermodynamic functions change within tolerance
615 : :
616 : : else
617 : : convth=1
618 : : end if
619 : :
620 : : ! Update F,E,S,C and eventually write them if converged
621 [ + + ]: 6 : if(convth==1)then
622 : 4 : part1=.true.
623 : 4 : write(msg,'(a,a,a)') ch10,&
624 : 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)'
625 : 4 : call wrtout(iout,msg,'COLL')
626 : 4 : call wrtout(thermal_unit,msg,'COLL')
627 : 4 : msg = ' # (A mol-c is the abbreviation of a mole-cell, that is, the'
628 : 4 : call wrtout(iout,msg,'COLL')
629 : 4 : call wrtout(thermal_unit,msg,'COLL')
630 : 4 : msg = ' # number of Avogadro times the atoms in a unit cell)'
631 : 4 : call wrtout(iout,msg,'COLL')
632 : 4 : call wrtout(thermal_unit,msg,'COLL')
633 : :
634 : : write(msg, '(a,a,a)' )&
635 : 4 : & ' harmonic_thermo : thermodynamic functions have converged',ch10,&
636 : 8 : & ' see main output file ...'
637 : 4 : call wrtout(std_out,msg,'COLL')
638 : : end if
639 : :
640 [ + + ]: 77 : do itemper=1,ntemper
641 : 70 : free0(itemper)=free(itemper)
642 : 70 : energy0(itemper)=energy(itemper)
643 : 70 : entropy0(itemper)=entropy(itemper)
644 : 70 : spheat0(itemper)=spheat(itemper)
645 : :
646 [ + + ]: 77 : if(convth==1)then
647 : 40 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
648 : 40 : write(msg,'(es11.3,5es15.7)') tmp+tol8,&
649 : 40 : & Ha_eV*e_Cb*Avogadro*free(itemper),&
650 : 40 : & Ha_eV*e_Cb*Avogadro*energy(itemper),&
651 : 40 : & Ha_eV*e_Cb*Avogadro*entropy(itemper),&
652 : 40 : & Ha_eV*e_Cb*Avogadro*spheat(itemper),&
653 : 80 : & wme(itemper)
654 : 40 : call wrtout(iout,msg,'COLL')
655 : 40 : call wrtout(thermal_unit,msg,'COLL')
656 : : end if
657 : : end do
658 : : end if
659 : :
660 [ + + ]: 17 : if(convth==1)exit
661 : : end do
662 : : end if
663 : :
664 [ + + ]: 14 : if(.not.part2)then
665 : : ! Atomic temperature factor calculation
666 [ + + ]: 65 : do iwchan=nwchan,1,-1
667 [ + + ]: 55 : if (wgijcnv(iwchan))cycle
668 : :
669 : : ! Calculate gij(k|w) and Bij(k)
670 : : ! Check if the density of states changed by more than dostol
671 : 25 : gijsum =zero
672 : 25 : wgijcnv(iwchan)=.true.
673 [ + + ]: 25 : if (ngrids>1) then
674 [ + + ]: 95 : do iatom=1,natom
675 [ + + ]: 545 : do ij=1,6
676 : 450 : gijerr=zero
677 [ + + ]: 217020 : do ichan=1,nchan2(iwchan)
678 : 216570 : gijsum = gijsum + gij(ij,iatom,ichan,iwchan)
679 : : gijerr=gijerr&
680 : : & +abs(ggij(ij,iatom,ichan,iwchan)/gnorm&
681 : 217020 : & - gij(ij,iatom,ichan,iwchan))
682 : : end do
683 [ - + ]: 525 : if(gijerr>anaddb_dtset%dostol) then
684 : 0 : wgijcnv(iwchan)=.false.
685 : 0 : exit
686 : : end if
687 : : end do
688 : : end do
689 : : else
690 : : gijerr=0.d0
691 : : end if
692 : :
693 : : ! gij(k|w) is updated
694 : :
695 [ + + ]: 10935 : do ichan=1,nchan2(iwchan)
696 [ + + ]: 72725 : do iatom=1,natom
697 [ + + ]: 443440 : do ij=1,6
698 : 432530 : gij(ij,iatom,ichan,iwchan)=ggij(ij,iatom,ichan,iwchan)/(gnorm/(3*natom))
699 : : end do
700 : : !if (iwchan==1) write (200+iatom,'(I6,6(E20.10,2x))') ichan, gij(1:6,iatom,ichan,iwchan)
701 : : end do
702 : : end do
703 : :
704 : : ! Write gijerr for each q sampling and w width
705 : :
706 : 25 : write(msg,'(a,a,i3,3i6,f10.5,f10.5)') ch10,&
707 : 25 : & ' iwchan,igqpt(i),gijsum, gij error= ',&
708 : 50 : & iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijsum,gijerr+tol10
709 : 25 : call wrtout(std_out,msg,'COLL')
710 : :
711 : : ! If the generalized DOS with a channel width is newly converged,
712 : : ! print it out and calculate Bij(k).
713 [ + - ]: 35 : if(wgijcnv(iwchan)) then
714 : :
715 [ + + ]: 25 : if (ngrids==1) then
716 : 5 : write(msg,'(a,i5,a)') ' gij with channel width= ',iwchan,':'
717 : : else
718 : 20 : write(msg,'(a,i5,a)') ' gij with channel width= ',iwchan,' newly converged'
719 : : end if
720 : 25 : call wrtout(iout,msg,'COLL')
721 : :
722 : 25 : write(msg,'(a,2i3,3i6,f10.5)')'iatom,iwchan,igqpt2(i),gij error= ',&
723 : 50 : & iatom,iwchan,igqpt2(1),igqpt2(2),igqpt2(3),gijerr+tol10
724 : 25 : call wrtout(iout,msg,'COLL')
725 : :
726 [ + + ]: 275 : do itemper=1,ntemper
727 : :
728 : : ! Put zeroes for Bij(k)
729 [ + + ]: 1450 : do iatom=1,natom
730 [ + + ]: 8650 : do ij=1,6
731 : 7200 : bbij(ij,iatom,itemper)=0._dp
732 : 8400 : vij(ij,iatom,itemper)=0._dp
733 : : end do
734 : : end do
735 : :
736 : 250 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
737 : : ! tmp in K
738 [ + - ]: 250 : if (tmp < tol6) cycle
739 : :
740 : 250 : dosinc=dble(iwchan)
741 : : !
742 [ + + ]: 109375 : do ichan=1,nchan2(iwchan)
743 : : !
744 : : !$wovert= \hbar*w / 2kT$, dimensionless
745 : 109100 : wovert=dosinc*(dble(ichan)-half)/Ha_cmm1/(two*kb_HaK*tmp)
746 : 109100 : expm2x=exp(-two*wovert)
747 [ + + ]: 727250 : do iatom=1,natom
748 : : ! factor contains 1 / (2 omega)
749 : : factor=Ha_cmm1/(two*dosinc*(dble(ichan)-half)) &
750 : 617900 : & *(one+expm2x)/(one-expm2x) /amu(Crystal%typat(iatom))/amu_emass
751 : :
752 : : ! this becomes * 0.5 * omega for the velocities
753 : : factorv=(half*dosinc*(dble(ichan)-half)/Ha_cmm1) &
754 : 617900 : & *(one+expm2x)/(one-expm2x) /amu(Crystal%typat(iatom))/amu_emass
755 : :
756 [ + + ]: 4434400 : do ij=1,6
757 : 3707400 : bbij(ij,iatom,itemper)=bbij(ij,iatom,itemper) + factor*gij(ij,iatom,ichan,iwchan)
758 : 4325300 : vij(ij,iatom,itemper)=vij(ij,iatom,itemper) + factorv*gij(ij,iatom,ichan,iwchan)
759 : : end do
760 : : end do
761 : :
762 : : end do
763 : :
764 : : end do
765 : :
766 : : ! B matrix is now in atomic unit in the Cartesian coordinates.
767 : : ! Check if Bij(k) changed within tolerance.
768 : 25 : convth=1
769 [ + + ]: 25 : if (ngrids>1) then
770 [ + + ]: 112 : do itemper=1,ntemper
771 [ + + ]: 309 : do iatom=1,natom
772 [ + + ]: 1467 : do ij=1,6
773 : 1261 : diffbb=bbij(ij,iatom,itemper)-bij(ij,iatom,itemper)
774 : : !if (diffbb > 1d-10 .and. diffbb/bij(ij,iatom,itemper) > thmtol) then
775 : : ! write(msg,'(a)' )' harmonic_thermo : Bij changes are larger than thmtol '
776 : : ! call wrtout(std_out,msg,'COLL')
777 : : ! convth=0
778 : : !end if
779 [ + + ]: 1261 : if (diffbb > 1d-10) then
780 [ + - ]: 409 : if (bij(ij,iatom,itemper) /= zero) then
781 [ + + ]: 409 : if (diffbb/bij(ij,iatom,itemper) > thmtol) then
782 : 11 : write(msg,'(a)' )' harmonic_thermo : Bij changes are larger than thmtol '
783 : 11 : call wrtout(std_out,msg,'COLL')
784 : : convth=0
785 : : end if
786 : : end if
787 : : end if
788 : 206 : if(convth==0)exit
789 : : end do
790 [ + + ]: 309 : if(convth==0)exit
791 : : end do
792 [ + + ]: 112 : if(convth==0)exit
793 : : end do
794 : : end if
795 : :
796 [ - + + + : 8700 : bij=bbij ! save for next iteration
+ + + + ]
797 : :
798 : : !Update Bij(k) and write them. B matrix printed in angstrom^2
799 : : !TODO : get rid of this version in the log and output file. Prefer
800 : : !external files
801 [ + + ]: 25 : if (convth==1) then
802 : : write(msg, '(a,a,a)' )&
803 : 14 : & ' B matrix elements as a function of T',ch10,&
804 : 28 : & ' Angstrom^2, cartesian coordinates'
805 : 14 : call wrtout(std_out,msg,'COLL')
806 : 14 : call wrtout(iout,msg,'COLL')
807 : :
808 [ + + ]: 154 : do itemper=1,ntemper
809 : : ! tmp in K
810 : 140 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
811 [ + + ]: 784 : do iatom=1,natom
812 : : write(iout,'(2i3,es11.3,6es12.4)')&
813 : 630 : & iwchan,iatom,tmp+tol10,&
814 : 630 : & Bohr_Ang**2*bij(1,iatom,itemper)+tol10,&
815 : 630 : & Bohr_Ang**2*bij(2,iatom,itemper)+tol10,&
816 : 630 : & Bohr_Ang**2*bij(3,iatom,itemper)+tol10,&
817 : 630 : & Bohr_Ang**2*bij(4,iatom,itemper)+tol10,&
818 : 630 : & Bohr_Ang**2*bij(5,iatom,itemper)+tol10,&
819 : 1400 : & Bohr_Ang**2*bij(6,iatom,itemper)+tol10
820 : : end do ! end loop over natom
821 : : end do ! end loop over ntemper
822 : :
823 : : ! Mean square velocity matrix printed in angstrom^2/picosec^2
824 : : write(msg, '(a,a,a)' )&
825 : 14 : & ' <vel^2> matrix elements as a function of T',ch10,&
826 : 28 : & ' Angstrom^2/(picosec)^2, cartesian coordinates'
827 : 14 : call wrtout(std_out,msg,'COLL')
828 : 14 : call wrtout(iout,msg,'COLL')
829 : :
830 [ + + ]: 154 : do itemper=1,ntemper
831 : : ! tmp in K
832 : 140 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
833 [ + + ]: 784 : do iatom=1,natom
834 [ + + ]: 4410 : vij(:,iatom,itemper)=Bohr_Ang**2*vij(:,iatom,itemper)/(Time_Sec*1.0e12)**2
835 : : ! The following check zeros out <v^2> if it is very small, in order to
836 : : ! avoid numerical noise being interpreted by the automatic tests as
837 : : ! something real. Note also that we compare it in
838 : : ! absolute value, that's because if any of the phonon frequencies are
839 : : ! computed as negative, <v^2> can take a negative value.
840 [ + + ]: 4410 : do icomp=1, 6
841 [ + + ]: 4410 : if (abs(vij(icomp,iatom,itemper)) < 1.0e-12) vij(icomp,iatom,itemper)=zero
842 : : end do
843 : : write(iout,'(2i3,es11.3,6es12.4)')&
844 : 630 : & iwchan,iatom,tmp+tol10,&
845 : 630 : & vij(1,iatom,itemper),&
846 : 630 : & vij(2,iatom,itemper),&
847 : 630 : & vij(3,iatom,itemper),&
848 : 630 : & vij(4,iatom,itemper),&
849 : 630 : & vij(5,iatom,itemper),&
850 : 1400 : & vij(6,iatom,itemper)
851 : : end do ! end loop over natom
852 : : end do ! end loop over ntemper
853 : : end if ! end check on convergence
854 : :
855 : :
856 : : ! keep this one !!!!!!!!!!!!!!!!!!
857 : : if (convth==1) then
858 : : write(msg, '(a,a,a)' )&
859 : 14 : & '# B matrix elements as a function of T, for each atom, and smallest omega channel width',ch10,&
860 : 28 : & '# Angstrom^2, cartesian coordinates'
861 : 14 : call wrtout(bij_unit,msg,'COLL')
862 [ + + ]: 77 : do iatom=1,natom
863 : 63 : write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
864 : 63 : call wrtout(bij_unit,msg,'COLL')
865 [ + + ]: 707 : do itemper=1,ntemper
866 : : ! tmp in K
867 : 630 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*dble(itemper-1)
868 : : write(msg,'(es11.3,6es12.4)')&
869 : 630 : & tmp,&
870 : 630 : & Bohr_Ang**2*bij(1,iatom,itemper),&
871 : 630 : & Bohr_Ang**2*bij(2,iatom,itemper),&
872 : 630 : & Bohr_Ang**2*bij(3,iatom,itemper),&
873 : 630 : & Bohr_Ang**2*bij(4,iatom,itemper),&
874 : 630 : & Bohr_Ang**2*bij(5,iatom,itemper),&
875 : 1260 : & Bohr_Ang**2*bij(6,iatom,itemper)
876 : 693 : call wrtout(bij_unit,msg,'COLL')
877 : : end do ! end loop over ntemper
878 : : end do ! end loop over natom
879 : :
880 : : ! Mean square velocity matrix printed in angstrom^2/picosec^2
881 : : write(msg, '(a,a,a)' )&
882 : 14 : & '# <vel^2> matrix elements as a function of T, for each atom, and smallest channel width',ch10,&
883 : 28 : & '# Angstrom^2/(picosec)^2, cartesian coordinates'
884 : 14 : call wrtout(vij_unit,msg,'COLL')
885 : :
886 [ + + ]: 77 : do iatom=1,natom
887 : 63 : write(msg, '(2a,i10)' ) ch10, '# for atom ', iatom
888 : 63 : call wrtout(vij_unit,msg,'COLL')
889 [ + + ]: 707 : do itemper=1,ntemper
890 : : ! tmp in K
891 : 630 : tmp=anaddb_dtset%tempermin+anaddb_dtset%temperinc*float(itemper-1)
892 [ + + ]: 4410 : vij(:,iatom,itemper)=Bohr_Ang**2*vij(:,iatom,itemper)/(Time_Sec*1.0e12)**2
893 : :
894 : : ! The following check zeros out <v^2> if it is very small, in order to
895 : : ! avoid numerical noise being interpreted by the automatic tests as
896 : : ! something real. Note also that we compare it in
897 : : ! absolute value, that's because if any of the phonon frequencies are
898 : : ! computed as negative, <v^2> can take a negative value.
899 [ + + ]: 4410 : do icomp=1, 6
900 [ + + ]: 4410 : if (abs(vij(icomp,iatom,itemper)) < 1.0e-12) vij(icomp,iatom,itemper)=zero
901 : : end do
902 : : write(vij_unit,'(es11.3,6es12.4)')&
903 : 630 : & tmp,&
904 : 630 : & vij(1,iatom,itemper),&
905 : 630 : & vij(2,iatom,itemper),&
906 : 630 : & vij(3,iatom,itemper),&
907 : 630 : & vij(4,iatom,itemper),&
908 : 630 : & vij(5,iatom,itemper),&
909 : 1323 : & vij(6,iatom,itemper)
910 : : end do ! end loop over ntemper
911 : : end do ! end loop over natom
912 : : end if ! end check on convergence
913 : :
914 : : if(convth==1)part2=.true.
915 : :
916 : : end if ! End of test on wgijcnv
917 : : end do ! End of loop over iwchan
918 : : end if ! End of part2
919 : :
920 [ + + ]: 14 : if(part1.and.part2)exit
921 : :
922 [ - + ]: 12 : ABI_FREE(indqpt1)
923 [ - + ]: 12 : ABI_FREE(qpt2)
924 [ - + ]: 12 : ABI_FREE(spqpt2)
925 [ - + ]: 12 : ABI_FREE(wtq)
926 [ - + ]: 30 : ABI_FREE(wtq2)
927 : :
928 : : end do ! End of the Loop on the q point grids
929 : :
930 [ - + ]: 4 : ABI_FREE(bbij)
931 [ - + ]: 4 : ABI_FREE(bij)
932 [ - + ]: 4 : ABI_FREE(energy)
933 [ - + ]: 4 : ABI_FREE(energy0)
934 [ - + ]: 4 : ABI_FREE(entropy)
935 [ - + ]: 4 : ABI_FREE(entropy0)
936 [ - + ]: 4 : ABI_FREE(free)
937 [ - + ]: 4 : ABI_FREE(free0)
938 [ - + ]: 4 : ABI_FREE(gdos)
939 [ - + ]: 4 : ABI_FREE(gg)
940 [ - + ]: 4 : ABI_FREE(gg_sum)
941 [ - + ]: 4 : ABI_FREE(gg_rest)
942 [ - + ]: 4 : ABI_FREE(ggij)
943 [ - + ]: 4 : ABI_FREE(gij)
944 [ - + ]: 4 : ABI_FREE(nchan2)
945 [ - + ]: 4 : ABI_FREE(spheat)
946 [ - + ]: 4 : ABI_FREE(spheat0)
947 [ - + ]: 4 : ABI_FREE(vij)
948 [ - + ]: 4 : ABI_FREE(wgcnv)
949 [ - + ]: 4 : ABI_FREE(wgijcnv)
950 [ + + ]: 4 : if(allocated(indqpt1)) then
951 : 2 : ABI_FREE(indqpt1)
952 : : end if
953 [ + + ]: 4 : if(allocated(qpt2)) then
954 : 2 : ABI_FREE(qpt2)
955 : : end if
956 [ + + ]: 4 : if(allocated(spqpt2)) then
957 : 2 : ABI_FREE(spqpt2)
958 : : end if
959 [ + + ]: 4 : if(allocated(wtq)) then
960 : 2 : ABI_FREE(wtq)
961 : : end if
962 [ + + ]: 4 : if(allocated(wtq2)) then
963 : 2 : ABI_FREE(wtq2)
964 : : end if
965 [ - + ]: 4 : ABI_FREE(gw)
966 [ - + ]: 4 : ABI_FREE(wme)
967 : :
968 [ - + ]: 4 : if(.not.part1)then
969 : : write(msg, '(a,a,a,a,a,a,a,a,a)' )&
970 : 0 : & 'No thermodynamical function is printed out :',ch10,&
971 : 0 : & 'the tolerance level that was asked ',ch10,&
972 : 0 : & 'has not been match with the grids specified.',ch10,&
973 : 0 : & 'Action: in the input file, increase the resolution',ch10,&
974 : 0 : & 'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
975 : 0 : ABI_ERROR(msg)
976 : : end if
977 : :
978 [ + + ]: 4 : if(.not.part2)then
979 : : write(msg,'(a,a,a,a,a,a,a,a,a)')&
980 : 2 : & 'No atomic factor tensor is printed out :',ch10,&
981 : 2 : & 'the tolerance level that was asked ',ch10,&
982 : 2 : & 'has not been match with the grids specified.',ch10,&
983 : 2 : & 'Action: in the input file, increase the resolution',ch10,&
984 : 4 : & 'of grids ng2qpt, or decrease the accuracy requirement thmtol.'
985 : 2 : ABI_WARNING(msg)
986 : : end if
987 : :
988 : 4 : close (thermal_unit)
989 : 4 : close (bij_unit)
990 : 4 : close (vij_unit)
991 : :
992 : : end subroutine harmonic_thermo
993 : : !!***
994 : :
995 : : end module m_harmonic_thermo
996 : : !!***
|
