Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 :
8 : module m_tdep_phi3
9 :
10 : use defs_basis
11 : use m_errors
12 : use m_abicore
13 : use m_numeric_tools
14 : use m_linalg_interfaces
15 : use m_xmpi
16 : use m_io_tools
17 : use m_crystal, only : crystal_t
18 : use m_ddb, only : ddb_type
19 : use m_ifc, only : ifc_type
20 : use m_tdep_abitypes, only : Qbz_type
21 : use m_htetra
22 : use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt
23 : use m_tdep_qpt, only : Qpoints_type
24 : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
25 : use m_tdep_latt, only : Lattice_type
26 : use m_tdep_shell, only : Shell_type
27 : use m_tdep_sym, only : Symmetries_type
28 : use m_tdep_phi2, only : Eigen_type,tdep_init_eigen2nd,tdep_destroy_eigen2nd
29 : use m_tdep_sampling, only : tdep_Sampling_type
30 : use m_tdep_solver, only : tdep_Solver_type
31 : use m_tdep_model, only : tdep_Model_type
32 : use m_tdep_constraints, only : Constraints_type
33 :
34 : implicit none
35 :
36 : public :: tdep_calc_ftot3
37 : public :: tdep_calc_phi3ref
38 : public :: tdep_write_phi3
39 : public :: tdep_build_phi3_333
40 : public :: tdep_calc_gruneisen
41 : public :: tdep_write_gruneisen
42 : public :: tdep_calc_alpha_gamma
43 : !FB public :: tdep_calc_lifetime1
44 : ! public :: tdep_calc_lifetime2
45 :
46 : contains
47 :
48 : !====================================================================================================
49 14 : subroutine tdep_calc_ftot3(Model,Invar,Shell3at,ucart,Sym)
50 :
51 : type(tdep_Model_type),intent(inout) :: Model
52 : type(atdep_dataset_type),intent(in) :: Invar
53 : type(Shell_type),intent(in) :: Shell3at
54 : type(Symmetries_type),intent(in) :: Sym
55 : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
56 :
57 : integer :: iatom,jatom,katom,isym,itrans,ishell,iatshell
58 : integer :: ii,jj,kk,istep
59 14 : double precision, allocatable :: Phi3_333(:,:,:)
60 14 : double precision, allocatable :: ucart_blas(:)
61 14 : double precision, allocatable :: ftot3(:,:)
62 :
63 574 : ABI_MALLOC(Phi3_333,(3,3,3)) ; Phi3_333(:,:,:)=0.d0
64 152750 : ABI_MALLOC(ftot3,(3*Invar%natom,Invar%my_nstep)); ftot3(:,:)=0.d0
65 2554 : do iatom=1,Invar%natom
66 10882 : do ishell=1,Shell3at%nshell
67 : ! Build the 3x3x3 IFC of an atom in this shell
68 8328 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
69 88976 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
70 79388 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
71 79388 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
72 79388 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
73 79388 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
74 79388 : call tdep_build_phi3_333(isym,Model%Phi3(:,:,:,ishell),Phi3_333,Sym,itrans)
75 : ! Calculation of the force components (third order)
76 1675476 : do istep=1,Invar%my_nstep
77 6430428 : do ii=1,3
78 20640880 : do jj=1,3
79 61922640 : do kk=1,3
80 : ftot3(3*(iatom-1)+ii,istep)=ftot3(3*(iatom-1)+ii,istep)+&
81 57159360 : & Phi3_333(ii,jj,kk)*ucart(jj,jatom,istep)*ucart(kk,katom,istep)
82 : end do !kk
83 : end do !jj
84 : end do !ii
85 : end do !istep
86 : end do !iatshell
87 : end do !ishell
88 : end do !iatom
89 14 : ABI_FREE(Phi3_333)
90 152694 : ftot3(:,:)=ftot3(:,:)/2.d0
91 :
92 7676 : ABI_MALLOC(ucart_blas,(3*Invar%natom)); ucart_blas(:)=0.d0
93 294 : do istep=1,Invar%my_nstep
94 152680 : ucart_blas(:)=0.d0
95 51080 : do jatom=1,Invar%natom
96 203480 : do jj=1,3
97 203200 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
98 : end do
99 : end do
100 : call DGEMM('T','N',1,1,3*Invar%natom,1./3.d0,ftot3(:,istep),3*Invar%natom,ucart_blas,&
101 280 : & 3*Invar%natom,0.d0,Model%Phi3UiUjUk(istep),3*Invar%natom)
102 : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=&
103 152694 : & Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)-ftot3(:,istep)
104 : end do
105 14 : ABI_FREE(ucart_blas)
106 14 : ABI_FREE(ftot3)
107 :
108 14 : end subroutine tdep_calc_ftot3
109 :
110 : !=====================================================================================================
111 14 : subroutine tdep_calc_phi3ref(Solver,Shell3at,Phi3_ref)
112 :
113 : type(tdep_Solver_type),intent(in) :: Solver
114 : type(Shell_type),intent(in) :: Shell3at
115 : double precision, intent(inout) :: Phi3_ref(3,3,3,Shell3at%nshell)
116 :
117 : integer :: ishell,ncoeff,ncoeff_prev
118 : integer :: ii,jj,kk,kappa
119 14 : double precision, allocatable :: Phi3_coeff(:)
120 :
121 220 : ABI_CALLOC(Phi3_coeff,(Solver%ncoeff3rd))
122 192 : Phi3_coeff(:) = Solver%theta(Solver%ncoeff1st+Solver%ncoeff2nd+1:Solver%ncoeff1st+Solver%ncoeff2nd+Solver%ncoeff3rd)
123 :
124 63 : do ishell=1,Shell3at%nshell
125 : ! Build the 3x3x3 IFC per shell
126 49 : ncoeff =Shell3at%ncoeff(ishell)
127 49 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)
128 49 : kappa=0
129 196 : do ii=1,3
130 637 : do jj=1,3
131 1911 : do kk=1,3
132 1323 : kappa=kappa+1
133 6570 : Phi3_ref(ii,jj,kk,ishell)=sum(Shell3at%proj(kappa,1:ncoeff,ishell)*Phi3_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
134 : end do
135 : end do
136 : end do
137 : ! Remove the rounding errors before writing (for non regression testing purposes)
138 210 : do ii=1,3
139 637 : do jj=1,3
140 1911 : do kk=1,3
141 1764 : if (abs(Phi3_ref(ii,jj,kk,ishell)).lt.tol8) Phi3_ref(ii,jj,kk,ishell)=zero
142 : end do
143 : end do
144 : end do
145 : end do
146 14 : ABI_FREE(Phi3_coeff)
147 :
148 14 : end subroutine tdep_calc_phi3ref
149 :
150 : !=====================================================================================================
151 14 : subroutine tdep_write_phi3(distance,Invar,Phi3_ref,Shell3at,Sym)
152 :
153 : type(atdep_dataset_type),intent(in) :: Invar
154 : type(Symmetries_type),intent(in) :: Sym
155 : type(Shell_type),intent(in) :: Shell3at
156 : double precision, intent(in) :: distance(Invar%natom,Invar%natom,4)
157 : double precision, intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
158 :
159 : integer :: ishell,isym,jatom,katom
160 : integer :: iatref,jatref,katref,iatshell,itrans
161 : integer :: ii,jj
162 : double precision :: tmp1,tmp2,tmp3
163 14 : double precision, allocatable :: Phi3_333(:,:,:)
164 :
165 14 : write(Invar%stdout,*) ' '
166 14 : write(Invar%stdout,*) '#############################################################################'
167 14 : write(Invar%stdout,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
168 14 : write(Invar%stdout,*) '#############################################################################'
169 :
170 : ! Write the IFCs in the data.out file (with others specifications:
171 : ! number of atoms in a shell, Trace...)
172 574 : ABI_MALLOC(Phi3_333,(3,3,3)) ; Phi3_333(:,:,:)=0.d0
173 63 : do ishell=1,Shell3at%nshell
174 49 : iatref=Shell3at%iatref(ishell)
175 63 : if (Shell3at%neighbours(iatref,ishell)%n_interactions.ne.0) then
176 49 : jatref=Shell3at%jatref(ishell)
177 49 : katref=Shell3at%katref(ishell)
178 49 : write(Invar%stdout,'(a,i4,a,i4,a)') ' ======== NEW SHELL (ishell=',ishell,&
179 98 : & '): There are',Shell3at%neighbours(iatref,ishell)%n_interactions,' atoms on this shell'
180 626 : do iatshell=1,Shell3at%neighbours(iatref,ishell)%n_interactions
181 577 : jatom=Shell3at%neighbours(iatref,ishell)%atomj_in_shell(iatshell)
182 577 : katom=Shell3at%neighbours(iatref,ishell)%atomk_in_shell(iatshell)
183 577 : isym =Shell3at%neighbours(iatref,ishell)%sym_in_shell(iatshell)
184 577 : itrans=Shell3at%neighbours(iatref,ishell)%transpose_in_shell(iatshell)
185 577 : call tdep_build_phi3_333(isym,Phi3_ref(:,:,:,ishell),Phi3_333,Sym,itrans)
186 577 : write(Invar%stdout,'(a,i4,a,i4)') ' For iatcell=',iatref,' ,with type=',mod(iatref-1,Invar%natom_unitcell)+1
187 577 : write(Invar%stdout,'(a,i4,a,i4)') ' For jatom =',jatom ,' ,with type=',mod(jatom -1,Invar%natom_unitcell)+1
188 577 : write(Invar%stdout,'(a,i4,a,i4)') ' For katom =',katom ,' ,with type=',mod(katom -1,Invar%natom_unitcell)+1
189 2308 : do jj=1,3
190 1731 : if (jj==1) write(Invar%stdout,'(a)') ' Phi3^{alpha beta x}='
191 1731 : if (jj==2) write(Invar%stdout,'(a)') ' Phi3^{alpha beta y}='
192 1731 : if (jj==3) write(Invar%stdout,'(a)') ' Phi3^{alpha beta z}='
193 7501 : do ii=1,3
194 5193 : if (abs(Phi3_333(ii,1,jj)).lt.5.d-7) then
195 1443 : tmp1=0.d0
196 : else
197 3750 : tmp1=Phi3_333(ii,1,jj)
198 : end if
199 5193 : if (abs(Phi3_333(ii,2,jj)).lt.5.d-7) then
200 1413 : tmp2=0.d0
201 : else
202 3780 : tmp2=Phi3_333(ii,2,jj)
203 : end if
204 5193 : if (abs(Phi3_333(ii,3,jj)).lt.5.d-7) then
205 1445 : tmp3=0.d0
206 : else
207 3748 : tmp3=Phi3_333(ii,3,jj)
208 : end if
209 6924 : write(Invar%stdout,'(2x,3(f9.6,1x))') tmp1,tmp2,tmp3
210 : end do
211 : end do
212 577 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (i,j) vector components:', (distance(iatref,jatom,jj+1),jj=1,3)
213 577 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (j,k) vector components:', (distance(jatom ,katom,jj+1),jj=1,3)
214 577 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (k,i) vector components:', (distance(katom,iatref,jj+1),jj=1,3)
215 626 : write(Invar%stdout,*) ' '
216 : end do !iatshell
217 : end if !n_interactions
218 : end do !ishell
219 14 : ABI_FREE(Phi3_333)
220 :
221 14 : end subroutine tdep_write_phi3
222 :
223 : !=====================================================================================================
224 10493 : subroutine tdep_calc_gruneisen(distance,Eigen2nd,Gruneisen,iqpt,Invar,Phi3_ref,qpt_cart,Rlatt_cart,Shell3at,Sym)
225 :
226 : type(Symmetries_type),intent(in) :: Sym
227 : type(atdep_dataset_type),intent(in) :: Invar
228 : type(Shell_type),intent(in) :: Shell3at
229 : type(Eigen_type),intent(in) :: Eigen2nd
230 : integer,intent(in) :: iqpt
231 : double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
232 : double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
233 : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
234 : double precision,intent(in) :: qpt_cart(3)
235 : double complex ,intent(out) :: Gruneisen(3*Invar%natom_unitcell,3,3)
236 :
237 : integer :: ii,jj,kk,iatcell,jatom,katom,itypat,jtypat,natom_unitcell,ll
238 : integer :: imode,nmode,jat_mod,jatcell,iatshell,ishell,isym,itrans
239 : double precision :: phase
240 : double precision :: F_strain(3,3)
241 10493 : double precision, allocatable :: Phi3_333(:,:,:)
242 10493 : double complex, allocatable :: eigen_prod(:,:,:,:,:),eigenvec_loc(:,:),omega(:)
243 :
244 : ! Define quantities
245 10493 : natom_unitcell=Invar%natom_unitcell
246 10493 : nmode=3*Invar%natom_unitcell
247 41972 : do ii=1,3
248 136409 : do jj=1,3
249 125916 : if (ii.eq.jj) then
250 : !FB F_strain(ii,jj)= 1./dsqrt(6.d0)
251 : !FB F_strain(ii,jj)=1./dsqrt(15.d0)
252 31479 : F_strain(ii,jj)= 1.d0
253 : !FB F_strain(ii,jj)= 0.8
254 : else
255 : !FB F_strain(ii,jj)= 1./dsqrt(6.d0)/2.d0
256 : !FB F_strain(ii,jj)=2./dsqrt(15.d0)
257 62958 : F_strain(ii,jj)=0.d0
258 : !FB F_strain(ii,jj)= 0.3
259 : end if
260 : end do
261 : end do
262 :
263 : ! Calculation of the eigenvectors product (and phonon frequencies)
264 4146685 : ABI_MALLOC(eigen_prod, (3,3,natom_unitcell,natom_unitcell,nmode)) ; eigen_prod(:,:,:,:,:)=czero
265 513241 : ABI_MALLOC(eigenvec_loc,(3*natom_unitcell,3*natom_unitcell)) ; eigenvec_loc(:,:) =czero
266 31275 : do iatcell=1,natom_unitcell
267 75545 : do jatcell=1,natom_unitcell
268 197862 : do ii=1,3
269 575510 : do jj=1,3
270 : eigenvec_loc((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
271 531240 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
272 : end do !ii
273 : end do !jj
274 : end do !jatcell
275 : end do !iatcell
276 93825 : ABI_MALLOC(omega,(3*natom_unitcell)); omega(:)=czero
277 31275 : do iatcell=1,natom_unitcell
278 75545 : do jatcell=1,natom_unitcell
279 197862 : do ii=1,3
280 575510 : do jj=1,3
281 3230592 : do imode=1,nmode
282 : eigen_prod(ii,jj,iatcell,jatcell,imode)=conjg(eigenvec_loc((iatcell-1)*3+ii,imode))*&
283 2699352 : & eigenvec_loc((jatcell-1)*3+jj,imode)
284 : omega(imode)=omega(imode)+eigen_prod(ii,jj,iatcell,jatcell,imode)*&
285 : & dcmplx(Eigen2nd%dynmat(1,ii,iatcell,jj,jatcell,iqpt),&
286 3097782 : & Eigen2nd%dynmat(2,ii,iatcell,jj,jatcell,iqpt))
287 : !FB write(Invar%stdlog,'(a,5(1x,i4),2(1x,f15.3))') 'imode,iatcell,jatcell,ii,jj=',imode,iatcell,jatcell,ii,jj
288 : !FB eigen_prod(ii,jj,iatcell,jatcell,imode)
289 : end do !imode
290 : end do !ii
291 : end do !jj
292 : end do !jatcell
293 : end do !iatcell
294 :
295 72839 : do imode=1,nmode
296 72839 : if (abs(aimag(omega(imode))).gt.tol8) then
297 0 : write(Invar%stdlog,'(a,1x,e15.8,1x,a,i4)') '>>>> WARNING : Imaginary part of the phonon frequency is not zero (',&
298 0 : & aimag(omega(imode)),') for mode :',imode
299 : !FB stop
300 : end if
301 : end do
302 : !FB write(Invar%stdlog,'(a,1(1x,i4),6(1x,f15.3))') 'iqpt,Grun=',iqpt,(dsqrt(real(omega(imode))),imode=1,6)
303 10493 : ABI_FREE(omega)
304 :
305 :
306 : ! Calculation of the Gruneisen
307 419720 : ABI_MALLOC(Phi3_333,(3,3,3)) ; Phi3_333(:,:,:)=0.d0
308 697523 : Gruneisen(:,:,:)=czero
309 42430 : do ishell=1,Shell3at%nshell
310 113375 : do iatcell=1,Invar%natom_unitcell
311 70945 : itypat=Invar%typat_unitcell(iatcell)
312 811660 : do iatshell=1,Shell3at%neighbours(iatcell,ishell)%n_interactions
313 708778 : jatom=Shell3at%neighbours(iatcell,ishell)%atomj_in_shell(iatshell)
314 708778 : katom=Shell3at%neighbours(iatcell,ishell)%atomk_in_shell(iatshell)
315 708778 : isym =Shell3at%neighbours(iatcell,ishell)%sym_in_shell(iatshell)
316 708778 : itrans=Shell3at%neighbours(iatcell,ishell)%transpose_in_shell(iatshell)
317 708778 : jat_mod=mod(jatom+Invar%natom_unitcell-1,Invar%natom_unitcell)+1
318 708778 : jtypat=Invar%typat_unitcell(jat_mod)
319 708778 : phase=0.d0
320 2835112 : do jj=1,3
321 2835112 : phase=phase+2*pi*Rlatt_cart(jj,iatcell,jatom)*qpt_cart(jj)
322 : end do
323 708778 : call tdep_build_phi3_333(isym,Phi3_ref(:,:,:,ishell),Phi3_333,Sym,itrans)
324 2906057 : do ii=1,3
325 9214114 : do jj=1,3
326 27642342 : do kk=1,3
327 82927026 : do ll=1,3
328 476827830 : do imode=1,nmode
329 : Gruneisen(imode,kk,ll)=Gruneisen(imode,kk,ll)-dcmplx(Phi3_333(ii,jj,kk),0.d0)*&
330 : & eigen_prod(ii,jj,iatcell,jat_mod,imode)*&
331 : & exp(dcmplx(0.d0,phase))*dcmplx(F_strain(kk,ll),0.d0)*dcmplx(distance(iatcell,katom,ll+1),0.d0)/&
332 : & dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)/&
333 457690824 : & dcmplx(6*Eigen2nd%eigenval(imode,iqpt)**2,0.d0)
334 : !FB do ll=1,3
335 : !FB Gruneisen(imode,kk,ll)=Gruneisen(imode,kk,ll)-dcmplx(Phi3_333(ii,jj,kk),0.d0)*&
336 : !FB& eigen_prod(ii,jj,iatcell,jat_mod,imode)*&
337 : !FB& exp(dcmplx(0.d0,phase))*F_strain(kk,ll)*dcmplx(distance(iatcell,katom,ll+1),0.d0)/&
338 : !FB& dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)/&
339 : !FB& dcmplx(6*Eigen2nd%eigenval(imode,iqpt)**2,0.d0)
340 : !FB!FB write(Invar%stdlog,'(a,4(1x,i4),2(1x,f15.3))') 'ii,jj,kk,imode,Grun=',ii,jj,kk,imode,Gruneisen(imode,kk,kk)
341 : !FB end do !ll
342 : end do !imode
343 : end do !ll
344 : end do !kk
345 : end do !jj
346 : end do !ii
347 : end do !iatshell
348 : end do !iatcell
349 : end do !ishell
350 :
351 : !FB Gruneisen(:,:,:)=Gruneisen(:,:,:)*dcmplx(dsqrt(6.d0),0.d0)
352 : !FB do imode=1,nmode
353 : !FB write(Invar%stdlog,'(a,1(1x,i4),6(1x,f15.3))') 'imode,Grun=',imode,(Gruneisen(imode,kk,kk),kk=1,3)
354 : !FB end do
355 10493 : ABI_FREE(Phi3_333)
356 10493 : ABI_FREE(eigen_prod)
357 10493 : ABI_FREE(eigenvec_loc)
358 :
359 10493 : end subroutine tdep_calc_gruneisen
360 :
361 : !=====================================================================================================
362 868131 : subroutine tdep_build_phi3_333(isym,Phi3_ref,Phi3_333,Sym,itrans)
363 :
364 : type(Symmetries_type),intent(in) :: Sym
365 : double precision, intent(in) :: Phi3_ref(3,3,3)
366 : double precision, intent(out) :: Phi3_333(3,3,3)
367 : integer,intent(in) :: isym,itrans
368 :
369 : integer :: alpha,beta,gama
370 : integer :: ii,jj,kk,ee,ff,gg,mu,nu,xi
371 : double precision :: Phi3_tmp(3,3,3)
372 :
373 :
374 : ! Transform in the new basis wrt S_ref
375 868131 : Phi3_333(:,:,:)=zero
376 3472524 : do mu=1,3
377 11285703 : do alpha=1,3
378 33857109 : do nu=1,3
379 101571327 : do beta=1,3
380 304713981 : do xi=1,3
381 914141943 : do gama=1,3
382 : Phi3_333(mu,nu,xi)=Phi3_333(mu,nu,xi)+&
383 843823332 : & Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu,beta,isym,1)*Sym%S_ref(xi,gama,isym,1)*Phi3_ref(alpha,beta,gama)
384 : end do
385 : end do
386 : end do
387 : end do
388 : end do
389 : end do
390 :
391 : ! Take into account the 6 allowed permutations
392 868131 : Phi3_tmp(:,:,:)=Phi3_333(:,:,:)
393 868131 : if ((itrans.lt.1).or.(itrans.gt.6)) then
394 0 : ABI_BUG('This value of the symmetry index is not permitted')
395 : end if
396 3472524 : do ii=1,3
397 11285703 : do jj=1,3
398 33857109 : do kk=1,3
399 23439537 : if (itrans==1) then ; ee=ii ; ff=jj ; gg=kk ; endif !\Phi3_efg
400 23439537 : if (itrans==2) then ; ee=ii ; ff=kk ; gg=jj ; endif !\Phi3_egf
401 23439537 : if (itrans==3) then ; ee=jj ; ff=ii ; gg=kk ; endif !\Phi3_feg
402 23439537 : if (itrans==4) then ; ee=jj ; ff=kk ; gg=ii ; endif !\Phi3_fge
403 23439537 : if (itrans==5) then ; ee=kk ; ff=ii ; gg=jj ; endif !\Phi3_gef
404 23439537 : if (itrans==6) then ; ee=kk ; ff=jj ; gg=ii ; endif !\Phi3_gfe
405 31252716 : Phi3_333(ee,ff,gg)=Phi3_tmp(ii,jj,kk)
406 : end do
407 : end do
408 : end do
409 :
410 868131 : end subroutine tdep_build_phi3_333
411 :
412 : !=====================================================================================================
413 14 : subroutine tdep_calc_alpha_gamma(distance,Eigen2nd,Invar,Lattice,MPIdata,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym)
414 :
415 : type(Eigen_type),intent(in) :: Eigen2nd
416 : type(atdep_dataset_type),intent(in) :: Invar
417 : type(Lattice_type),intent(inout) :: Lattice
418 : type(Shell_type),intent(in) :: Shell3at
419 : type(Symmetries_type),intent(in) :: Sym
420 : type(Qbz_type),intent(in) :: Qbz
421 : type(MPI_enreg_type), intent(in) :: MPIdata
422 : double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
423 : double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
424 : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
425 :
426 : integer :: iq_ibz,ii,jj,nmode,ntemp,itemp,kk,imode,ll
427 : integer :: alpha,beta
428 : double precision :: k_B,wovert,xx,C_v,E_th,Vp,Vs,P_th2
429 14 : double precision, allocatable :: heatcapa(:),grun_thermo(:,:,:),u_vib(:)
430 14 : double precision, allocatable :: heatcapa_HA(:,:),grun_thermo_HA(:,:,:,:),p_thermo1(:,:,:)
431 14 : double precision, allocatable :: p_thermo_HA(:,:,:,:),u_vib_HA(:,:)
432 14 : double precision, allocatable :: p_thermo2(:),tmp(:)
433 : double precision :: Gama_tensor(3,3),alpha_v_tensor(3,3),P_th_tensor(3,3),P_th1_tensor(3,3)
434 14 : double complex, allocatable :: Gruneisen(:,:,:),Grun_mean(:)
435 : character(len=500) :: message
436 : integer :: ierr
437 :
438 14 : ierr = 0
439 :
440 1008 : ABI_MALLOC(Gruneisen,(3*Invar%natom_unitcell,3,3)); Gruneisen(:,:,:)=czero
441 126 : ABI_MALLOC(Grun_mean,(3*Invar%natom_unitcell)) ; Grun_mean(:) =czero
442 :
443 14 : ntemp=1000
444 14 : nmode=3*Invar%natom_unitcell
445 14 : k_B=kb_HaK*Ha_eV
446 126 : ABI_MALLOC(heatcapa,(nmode)) ; heatcapa (:) =0.d0
447 980 : ABI_MALLOC(grun_thermo,(nmode,3,3)); grun_thermo(:,:,:)=0.d0
448 952 : ABI_MALLOC(p_thermo1,(nmode,3,3)) ; p_thermo1 (:,:,:)=0.d0
449 14014 : ABI_MALLOC(p_thermo2,(ntemp)) ; p_thermo2 (:) =0.d0
450 112 : ABI_MALLOC(u_vib,(nmode)) ; u_vib (:) =0.d0
451 98056 : ABI_MALLOC(heatcapa_HA,(nmode,ntemp)) ; heatcapa_HA (:,:) =0.d0
452 882252 : ABI_MALLOC(grun_thermo_HA,(nmode,ntemp,3,3)); grun_thermo_HA(:,:,:,:)=0.d0
453 882210 : ABI_MALLOC(p_thermo_HA,(nmode,ntemp,3,3)) ; p_thermo_HA (:,:,:,:)=0.d0
454 98042 : ABI_MALLOC(u_vib_HA,(nmode,ntemp)) ; u_vib_HA (:,:) =0.d0
455 : ! Loop over irreducible q-points
456 : ! =======================
457 1978 : do iq_ibz=1,Qbz%nqibz
458 :
459 : ! Compute the gruneisen for this q-point
460 : ! ======================================
461 7856 : if ((sum(abs(Qbz%qibz_cart(:,iq_ibz)))).lt.tol8) cycle ! G point
462 : call tdep_calc_gruneisen(distance,Eigen2nd,Gruneisen,iq_ibz,Invar,&
463 1962 : & Phi3_ref,Qbz%qibz_cart(:,iq_ibz),Rlatt_cart,Shell3at,Sym)
464 13086 : Grun_mean(:)=czero
465 13086 : do ii=1,3*Invar%natom_unitcell
466 46458 : do jj=1,3
467 144612 : do kk=1,3
468 133488 : Grun_mean(ii)=Grun_mean(ii)+Gruneisen(ii,jj,kk)
469 : end do
470 : end do
471 : end do
472 13086 : if (sum(abs(aimag(Grun_mean(:)))).gt.tol8) then
473 0 : write(message,'(i5,1x,100(e15.6,1x))') iq_ibz,( real(Grun_mean(ii)),ii=1,nmode)
474 0 : ABI_ERROR_NOSTOP(message,ierr)
475 0 : write(message,'(i5,1x,100(e15.6,1x))') iq_ibz,(aimag(Grun_mean(ii)),ii=1,nmode)
476 0 : ABI_ERROR_NOSTOP(message,ierr)
477 0 : ABI_ERROR('tdep_calc_alpha_gamma : The imaginary part of the Gruneisen is not equal to zero')
478 : end if
479 : ! If the Gruneisen is isotropic
480 : !FB do ii=1,3
481 : !FB tmp1=0.d0
482 : !FB do jj=1,3
483 : !FB tmp1=tmp1 + real(Gruneisen(ii,jj,jj))
484 : !FB end do
485 : !FB do jj=1,3
486 : !FB Gruneisen(ii,jj,jj)=dcmplx(tmp1/3.d0,0.d0)
487 : !FB end do
488 : !FB write(Invar%stdlog,'(2(i5,1x),100(e15.6,1x))') iq_ibz,ii,Gruneisen(ii,jj,jj)
489 : !FB end do
490 :
491 : ! Compute the heat capacity and thermodynamical gruneisen parameter at present temperature
492 : ! ========================================================================================
493 1962 : wovert=1.d0/(2*Invar%temperature*k_B)
494 13086 : do ii=1,nmode
495 11124 : xx=Eigen2nd%eigenval(ii,iq_ibz)*Ha_eV
496 11124 : if (xx.le.0) cycle
497 11124 : heatcapa(ii) =heatcapa(ii) +(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)
498 : ! Gamma= sum_i Gamma_i*C_Vi / C_V
499 144612 : grun_thermo(ii,:,:)=grun_thermo(ii,:,:)+(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)*real(Gruneisen(ii,:,:))*3.d0
500 : ! P = sum_i Gamma_i*E_i / V
501 144612 : p_thermo1(ii,:,:) =p_thermo1(ii,:,:) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)*real(Gruneisen(ii,:,:))*3.d0
502 13086 : u_vib(ii) =u_vib(ii) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)
503 : end do
504 : ! Compute the heat capacity and thermodynamical gruneisen parameter as a function of temperature
505 : ! ==============================================================================================
506 1963976 : do itemp=1,ntemp
507 1962000 : wovert=1.d0/(2*real(itemp)*10*k_B)
508 13087964 : do ii=1,nmode
509 11124000 : xx=Eigen2nd%eigenval(ii,iq_ibz)*Ha_eV
510 11124000 : if (xx.le.0) cycle
511 11124000 : heatcapa_HA(ii,itemp) =heatcapa_HA(ii,itemp) +(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)
512 : grun_thermo_HA(ii,itemp,:,:)=grun_thermo_HA(ii,itemp,:,:)+(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)&
513 144612000 : & *real(Gruneisen(ii,:,:))*3.d0
514 : p_thermo_HA(ii,itemp,:,:) =p_thermo_HA(ii,itemp,:,:) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)&
515 144612000 : & *real(Gruneisen(ii,:,:))*3.d0
516 13086000 : u_vib_HA(ii,itemp) =u_vib_HA(ii,itemp) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)
517 : end do
518 : end do
519 : end do
520 14 : ABI_FREE(Gruneisen)
521 14 : ABI_FREE(Grun_mean)
522 : ! Compute the pressure as the integral of Gamma*C_v overt T
523 14014 : ABI_MALLOC(tmp,(ntemp)) ; tmp(:)=0.d0
524 14014 : do itemp=1,ntemp
525 56014 : do ii=1,3
526 182000 : do jj=1,3
527 924000 : tmp(itemp)=tmp(itemp)+sum(grun_thermo_HA(:,itemp,ii,jj))/3.d0
528 : end do
529 : end do
530 : end do
531 14 : call simpson_int(ntemp,10.d0,tmp,p_thermo2)
532 14 : ABI_FREE(tmp)
533 :
534 14 : if (MPIdata%iam_master) then
535 14 : open(unit=20,file=trim(Invar%output_prefix)//'_thermo3.dat')
536 14 : open(unit=21,file=trim(Invar%output_prefix)//'_alpha_gamma.dat')
537 14 : write(20,'(a)')'# T(K) C_v(k_B/fu) Gamma alpha_v*10^6(K^-1) E_th(eV) P_th(GPa)'
538 14 : write(20,'(a,72x,a)')'#',' ----------------------------------------------'
539 14 : write(20,'(a,72x,a)')'#',' {sum G_i.U_iV} {int G.C_v/V dT} {G.U/V}'
540 14 : write(21,'(a)')'# T(K) Gamma_11 Gamma_22 Gamma_33 alpha_11 alpha_22 alpha_33',&
541 28 : & ' alpha_12 alpha_13 alpha_23'
542 14 : write(20,'(2a)')'# --------------------------------------------------------------------------------------------------',&
543 28 : '--------------------'
544 : end if
545 14014 : do itemp=1,ntemp
546 98000 : C_v =sum(heatcapa_HA (:,itemp))
547 98000 : E_th =sum(u_vib_HA (:,itemp))
548 14000 : alpha_v_tensor(:,:)=zero
549 14000 : Gama_tensor (:,:)=zero
550 14000 : P_th1_tensor (:,:)=zero
551 14000 : P_th_tensor (:,:)=zero
552 98000 : do imode=1,nmode
553 1092000 : Gama_tensor (:,:)=Gama_tensor (:,:)+grun_thermo_HA(imode,itemp,:,:)/C_v
554 1106000 : P_th1_tensor(:,:)=P_th1_tensor(:,:)+p_thermo_HA (imode,itemp,:,:)/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
555 : end do
556 : ! alpha=sum{I=1,6;J=1,6} S(I,J) Gama_tensor (J)
557 56000 : do ii=1,3
558 182000 : do jj=1,3
559 126000 : if (ii.eq.1.and.jj.eq.1) alpha=1
560 126000 : if (ii.eq.2.and.jj.eq.2) alpha=2
561 126000 : if (ii.eq.3.and.jj.eq.3) alpha=3
562 126000 : if (ii.eq.2.and.jj.eq.3) alpha=4
563 126000 : if (ii.eq.1.and.jj.eq.3) alpha=5
564 126000 : if (ii.eq.1.and.jj.eq.2) alpha=6
565 126000 : if (ii.eq.3.and.jj.eq.2) alpha=4
566 126000 : if (ii.eq.3.and.jj.eq.1) alpha=5
567 126000 : if (ii.eq.2.and.jj.eq.1) alpha=6
568 504000 : do kk=1,3
569 1638000 : do ll=1,3
570 1134000 : if (kk.eq.1.and.ll.eq.1) beta=1
571 1134000 : if (kk.eq.2.and.ll.eq.2) beta=2
572 1134000 : if (kk.eq.3.and.ll.eq.3) beta=3
573 1134000 : if (kk.eq.2.and.ll.eq.3) beta=4
574 1134000 : if (kk.eq.1.and.ll.eq.3) beta=5
575 1134000 : if (kk.eq.1.and.ll.eq.2) beta=6
576 1134000 : if (kk.eq.3.and.ll.eq.2) beta=4
577 1134000 : if (kk.eq.3.and.ll.eq.1) beta=5
578 1134000 : if (kk.eq.2.and.ll.eq.1) beta=6
579 : alpha_v_tensor(ii,jj)=alpha_v_tensor(ii,jj)+Lattice%Sij(alpha,beta)*Gama_tensor(kk,ll)*C_v*kb_HaK*Ha_J/1.d9/&
580 1512000 : & (Lattice%ucvol*Bohr_Ang**3*1.d-30)
581 : end do
582 : end do
583 168000 : P_th_tensor(ii,jj) =Gama_tensor(ii,jj)*E_th/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
584 : end do
585 : end do
586 14000 : P_th2 = p_thermo2(itemp)*k_B/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
587 14014 : if (MPIdata%iam_master) then
588 14000 : write(20,'(1x,i5,7(1x,f15.3))') itemp*10,C_v,(Gama_tensor(1,1)+Gama_tensor(2,2)+Gama_tensor(3,3))/3.d0,&
589 14000 : & (alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*1.d6,&
590 14000 : & E_th,(P_th1_tensor(1,1)+P_th1_tensor(2,2)+P_th1_tensor(3,3))/3.d0,&
591 28000 : & P_th2,(P_th_tensor(1,1)+P_th_tensor(2,2)+P_th_tensor(3,3))/3.d0
592 14000 : write(21,'(1x,i5,9(1x,f15.3))') itemp*10,Gama_tensor(1,1),Gama_tensor(2,2),Gama_tensor(3,3),&
593 14000 : & alpha_v_tensor(1,1)*1.d6,alpha_v_tensor(2,2)*1.d6,alpha_v_tensor(3,3)*1.d6,&
594 28000 : & alpha_v_tensor(1,2)*1.d6,alpha_v_tensor(1,3)*1.d6,alpha_v_tensor(2,3)*1.d6
595 : end if
596 : end do
597 14 : if (MPIdata%iam_master) then
598 14 : close(20)
599 14 : close(21)
600 : end if
601 14 : ABI_FREE(heatcapa_HA)
602 14 : ABI_FREE(grun_thermo_HA)
603 14 : ABI_FREE(p_thermo_HA)
604 14 : ABI_FREE(u_vib_HA)
605 :
606 : ! Summary
607 : ! =======
608 14 : write(Invar%stdout,*) ' '
609 14 : write(Invar%stdout,*) '#############################################################################'
610 14 : write(Invar%stdout,*) '####### Gruneisen parameter, Thermal expansion, Thermal pressure... #########'
611 14 : write(Invar%stdout,*) '#############################################################################'
612 98 : C_v =sum(heatcapa(:))
613 98 : E_th =sum(u_vib (:))
614 14 : alpha_v_tensor(:,:)=zero
615 14 : Gama_tensor (:,:)=zero
616 14 : P_th1_tensor (:,:)=zero
617 14 : P_th_tensor (:,:)=zero
618 98 : do imode=1,nmode
619 1092 : Gama_tensor (:,:)=Gama_tensor (:,:)+grun_thermo(imode,:,:)/C_v
620 1106 : P_th1_tensor(:,:)=P_th1_tensor(:,:)+p_thermo1 (imode,:,:)/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
621 : end do
622 : ! alpha=sum{I=1,6;J=1,6} S(I,J) Gama_tensor (J)
623 56 : do ii=1,3
624 182 : do jj=1,3
625 126 : if (ii.eq.1.and.jj.eq.1) alpha=1
626 126 : if (ii.eq.2.and.jj.eq.2) alpha=2
627 126 : if (ii.eq.3.and.jj.eq.3) alpha=3
628 126 : if (ii.eq.2.and.jj.eq.3) alpha=4
629 126 : if (ii.eq.1.and.jj.eq.3) alpha=5
630 126 : if (ii.eq.1.and.jj.eq.2) alpha=6
631 126 : if (ii.eq.3.and.jj.eq.2) alpha=4
632 126 : if (ii.eq.3.and.jj.eq.1) alpha=5
633 126 : if (ii.eq.2.and.jj.eq.1) alpha=6
634 504 : do kk=1,3
635 1638 : do ll=1,3
636 1134 : if (kk.eq.1.and.ll.eq.1) beta=1
637 1134 : if (kk.eq.2.and.ll.eq.2) beta=2
638 1134 : if (kk.eq.3.and.ll.eq.3) beta=3
639 1134 : if (kk.eq.2.and.ll.eq.3) beta=4
640 1134 : if (kk.eq.1.and.ll.eq.3) beta=5
641 1134 : if (kk.eq.1.and.ll.eq.2) beta=6
642 1134 : if (kk.eq.3.and.ll.eq.2) beta=4
643 1134 : if (kk.eq.3.and.ll.eq.1) beta=5
644 1134 : if (kk.eq.2.and.ll.eq.1) beta=6
645 : alpha_v_tensor(ii,jj)=alpha_v_tensor(ii,jj)+Lattice%Sij(alpha,beta)*Gama_tensor(kk,ll)*C_v*kb_HaK*Ha_J/1.d9/&
646 1512 : & (Lattice%ucvol*Bohr_Ang**3*1.d-30)
647 : end do
648 : end do
649 168 : P_th_tensor(ii,jj)=Gama_tensor(ii,jj)*E_th/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
650 : end do
651 : end do
652 14 : P_th2 =p_thermo2(int(Invar%temperature/10))*k_B/(Lattice%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
653 : Lattice%BulkModulus_S=Lattice%BulkModulus_T*(1.+(alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*&
654 14 : & (Gama_tensor(1,1)+ Gama_tensor(2,2)+ Gama_tensor(3,3))/3.d0*Invar%temperature)
655 14 : Lattice%HeatCapa_P=C_v*Lattice%BulkModulus_S/Lattice%BulkModulus_T
656 14 : Vp=dsqrt(1.d9*(Lattice%BulkModulus_S+4.d0*Lattice%Shear/3.d0)/Lattice%Density)
657 14 : Vs=dsqrt(1.d9*Lattice%Shear/Lattice%Density)
658 14 : write(Invar%stdout,'(a)') ' See the gruneisen.dat, alpha_gamma.dat and thermo3.dat files'
659 14 : write(Invar%stdout,'(a)') ' '
660 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Gruneisen parameter : Gamma=',&
661 28 : & (Gama_tensor(1,1)+Gama_tensor(2,2)+Gama_tensor(3,3))/3.d0
662 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' The Gruneisen matrix is |',&
663 28 : & Gama_tensor(1,1),Gama_tensor(1,2),Gama_tensor(1,3),'|'
664 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
665 28 : & Gama_tensor(2,1),Gama_tensor(2,2),Gama_tensor(2,3),'|'
666 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
667 28 : & Gama_tensor(3,1),Gama_tensor(3,2),Gama_tensor(3,3),'|'
668 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Thermal expansion (K^{-1}*10^6) : alpha_v=',&
669 28 : & (alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*1.d6
670 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' The thermal expansion matrix is |',&
671 28 : & alpha_v_tensor(1,1)*1.d6,alpha_v_tensor(1,2)*1.d6,alpha_v_tensor(1,3)*1.d6,'|'
672 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
673 28 : & alpha_v_tensor(2,1)*1.d6,alpha_v_tensor(2,2)*1.d6,alpha_v_tensor(2,3)*1.d6,'|'
674 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
675 28 : & alpha_v_tensor(3,1)*1.d6,alpha_v_tensor(3,2)*1.d6,alpha_v_tensor(3,3)*1.d6,'|'
676 14 : write(Invar%stdout,'(a)') ' Thermal pressure (in GPa) : '
677 14 : write(Invar%stdout,'(a)') ' ------- w intrinsic effects and w ZPE --------'
678 14 : write(Invar%stdout,'(a,1x,f15.3)') ' P_th=sum_i Gamma_i*E_i/V =',&
679 28 : & (P_th1_tensor(1,1)+P_th1_tensor(2,2)+P_th1_tensor(3,3))/3.d0
680 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' The thermal pressure matrix is|',&
681 28 : & P_th1_tensor(1,1),P_th1_tensor(1,2),P_th1_tensor(1,3),'|'
682 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
683 28 : & P_th1_tensor(2,1),P_th1_tensor(2,2),P_th1_tensor(2,3),'|'
684 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
685 28 : & P_th1_tensor(3,1),P_th1_tensor(3,2),P_th1_tensor(3,3),'|'
686 14 : write(Invar%stdout,'(a)') ' ------- w intrinsic effects and w/o ZPE --------'
687 14 : write(Invar%stdout,'(a,1x,f15.3)') ' P_th=integ{Gamma*C_v/V dT} =',P_th2
688 14 : write(Invar%stdout,'(a)') ' ------- w/o intrinsic effects and w/o ZPE --------'
689 14 : write(Invar%stdout,'(a,1x,f15.3)') ' P_th=Gamma*E_th/V =',&
690 28 : & (P_th_tensor(1,1)+P_th_tensor(2,2)+P_th_tensor(3,3))/3.d0
691 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' The thermal pressure matrix is|',&
692 28 : & P_th_tensor(1,1),P_th_tensor(1,2),P_th_tensor(1,3),'|'
693 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
694 28 : & P_th_tensor(2,1),P_th_tensor(2,2),P_th_tensor(2,3),'|'
695 14 : write(Invar%stdout,'(a,3(1x,f8.3),a)') ' |',&
696 28 : & P_th_tensor(3,1),P_th_tensor(3,2),P_th_tensor(3,3),'|'
697 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Volume (bohr^3 per unit cell): V=',Lattice%ucvol
698 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Thermal energy (eV) : E_th=',E_th
699 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Heat capacity at constant V (k_B/f.u.): C_v=',C_v
700 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Heat capacity at constant P (k_B/f.u.): C_p=',Lattice%HeatCapa_P
701 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Isothermal Bulk Modulus (GPa): B_T=',Lattice%BulkModulus_T
702 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Isentropic Bulk Modulus (GPa): B_S=',Lattice%BulkModulus_S
703 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Longitudinal sound velocity (m.s-1): Vp=',Vp
704 14 : write(Invar%stdout,'(a,1x,f15.3)') ' Transverse sound velocity (m.s-1): Vs=',Vs
705 14 : ABI_FREE(heatcapa)
706 14 : ABI_FREE(grun_thermo)
707 14 : ABI_FREE(p_thermo1)
708 14 : ABI_FREE(p_thermo2)
709 14 : ABI_FREE(u_vib)
710 :
711 14 : end subroutine tdep_calc_alpha_gamma
712 :
713 : !=====================================================================================================
714 14 : subroutine tdep_write_gruneisen(distance,Eigen2nd,Invar,Phi3_ref,Qpt,Rlatt_cart,Shell3at,Sym)
715 :
716 : type(Symmetries_type),intent(in) :: Sym
717 : type(atdep_dataset_type),intent(in) :: Invar
718 : type(Shell_type),intent(in) :: Shell3at
719 : type(Eigen_type),intent(in) :: Eigen2nd
720 : type(Qpoints_type),intent(in) :: Qpt
721 : double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
722 : double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
723 : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
724 :
725 : integer :: iqpt,nmode,ii,jj,kk
726 : double precision :: qpt_cart(3)
727 14 : double complex, allocatable :: Gruneisen(:,:,:)
728 14 : double complex, allocatable :: Grun_mean(:)
729 : character(len=500) :: message
730 : integer :: ierr
731 :
732 14 : ierr = 0
733 :
734 14 : nmode=3*Invar%natom_unitcell
735 1008 : ABI_MALLOC(Gruneisen,(3*Invar%natom_unitcell,3,3)); Gruneisen(:,:,:)=czero
736 126 : ABI_MALLOC(Grun_mean,(3*Invar%natom_unitcell)) ; Grun_mean(:) =czero
737 14 : open(unit=53,file=trim(Invar%output_prefix)//'_gruneisen.dat')
738 14 : open(unit=54,file=trim(Invar%output_prefix)//'_gruneisen-ij.dat')
739 :
740 8573 : do iqpt=1,Qpt%nqpt
741 :
742 : ! Compute the Gruneisen
743 : ! =====================
744 59949 : Grun_mean(:)=czero
745 34236 : qpt_cart(:)=Qpt%qpt_cart(:,iqpt)
746 34236 : if ((sum(abs(Qpt%qpt_red(:,iqpt)))).lt.tol8) cycle ! G point
747 34132 : if (abs(sum(Qpt%qpt_red(:,iqpt)**2)-1.d0).lt.tol8) cycle ! Gp point
748 8531 : call tdep_calc_gruneisen(distance,Eigen2nd,Gruneisen,iqpt,Invar,Phi3_ref,qpt_cart,Rlatt_cart,Shell3at,Sym)
749 59753 : do ii=1,3*Invar%natom_unitcell
750 213419 : do jj=1,3
751 665886 : do kk=1,3
752 460998 : write(54,'(4(i5,1x),500(e15.6,1x))') iqpt,ii,jj,kk,Gruneisen(ii,jj,kk)
753 614664 : if (abs(aimag(Gruneisen(ii,jj,kk))).gt.tol8) then
754 0 : ABI_WARNING('Real part :')
755 0 : write(message,'(4(i5,1x),100(e15.6,1x))') iqpt,ii,jj,kk,real(Gruneisen(ii,jj,kk))
756 0 : ABI_WARNING(message)
757 0 : ABI_WARNING('Imaginary part :')
758 0 : write(message,'(4(i5,1x),100(e15.6,1x))') iqpt,ii,jj,kk,aimag(Gruneisen(ii,jj,kk))
759 0 : ABI_WARNING(message)
760 0 : ABI_WARNING('tdep_write_gruneisen : The imaginary part of the Gruneisen is not equal to zero')
761 : end if
762 : end do
763 : end do
764 : end do
765 59753 : do ii=1,3*Invar%natom_unitcell
766 213419 : do jj=1,3
767 665886 : do kk=1,3
768 614664 : Grun_mean(ii)=Grun_mean(ii)+Gruneisen(ii,jj,kk)
769 : end do
770 : end do
771 : end do
772 :
773 : ! Write the Gruneisen
774 : ! ===================
775 59767 : if (sum(abs(aimag(Grun_mean(:)))).gt.3*Invar%natom_unitcell*tol8) then
776 0 : ABI_WARNING('Real part :')
777 0 : write(message,'(i5,1x,100(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode)
778 0 : ABI_WARNING(message)
779 0 : ABI_WARNING('Imaginary part :')
780 0 : write(message,'(i5,1x,100(e15.6,1x))') iqpt,(aimag(Grun_mean(ii)),ii=1,nmode)
781 0 : ABI_WARNING(message)
782 0 : ABI_WARNING('tdep_write_gruneisen : The imaginary part of Grun_mean is not equal to zero')
783 : else
784 : !FB write(53,'(i5,1x,500(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode),&
785 : !FB ((real(Grun_shell(ii,jj)),ii=1,nmode),jj=1,Shell3at%nshell)
786 59753 : write(53,'(i5,1x,500(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode)
787 : end if
788 : end do
789 14 : close(53)
790 14 : close(54)
791 14 : ABI_FREE(Grun_mean)
792 14 : ABI_FREE(Gruneisen)
793 :
794 14 : end subroutine tdep_write_gruneisen
795 :
796 : !=====================================================================================================
797 : end module m_tdep_phi3
798 :
799 : !FB!=====================================================================================================
800 : !FBsubroutine tdep_calc_lifetime1(Crystal,distance,Eigen2nd,Ifc,Invar,Lattice,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym)
801 : !FB
802 : !FB type(crystal_t),intent(in) :: Crystal
803 : !FB type(Symmetries_type),intent(in) :: Sym
804 : !FB type(atdep_dataset_type),intent(in) :: Invar
805 : !FB type(Shell_type),intent(in) :: Shell3at
806 : !FB type(Lattice_type),intent(in) :: Lattice
807 : !FB type(Eigen_type),intent(in) :: Eigen2nd
808 : !FB type(Qbz_type),intent(in) :: Qbz
809 : !FB type(ifc_type),intent(in) :: Ifc
810 : !FB double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
811 : !FB double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
812 : !FB double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
813 : !FB
814 : !FB integer :: iq_bz,jq_bz,kq_bz,nmode,ii,tot_count
815 : !FB integer :: okp_count,okm_count
816 : !FB double precision :: iqbz(3),jqbz(3),kqbz(3),ipjpk(3),imjpk(3)
817 : !FB! double complex, allocatable :: Gruneisen(:,:,:)
818 : !FB! double complex, allocatable :: Grun_mean(:)
819 : !FB integer :: ierr
820 : !FB
821 : !FB ierr = 0
822 : !FB
823 : !FB nmode=3*Invar%natom_unitcell
824 : !FB tot_count=zero
825 : !FB okp_count=zero
826 : !FB okm_count=zero
827 : !FB do ii=1,3
828 : !FB write(Invar%stdout,'(a,3(e15.6,1x),1x)') 'Gmet=',Lattice%gmet(ii,:)
829 : !FB write(Invar%stdout,'(a,3(e15.6,1x),1x)') 'Gprim=',Lattice%gprim(ii,:)
830 : !FB end do
831 : !FB do iq_bz=1,Qbz%nqbz
832 : !FB iqbz(:)=Qbz%qbz(:,iq_bz)
833 : !FB! iqbz(:)=0.2
834 : !FB write(Invar%stdout,'(a,3(e15.6,1x),1x)') 'Qbz1=',iqbz(:)
835 : !FB do jq_bz=1,Qbz%nqbz
836 : !FB jqbz(:)=Qbz%qbz(:,jq_bz)
837 : !FB! jqbz(:)=0.3
838 : !FB do kq_bz=1,Qbz%nqbz
839 : !FB tot_count=tot_count+1
840 : !FB kqbz(:)=Qbz%qbz(:,kq_bz)
841 : !FB! kqbz(:)=0.5
842 : !FB ipjpk(:)=iqbz(:)+jqbz(:)+kqbz(:)
843 : !FB imjpk(:)=iqbz(:)-jqbz(:)+kqbz(:)
844 : !FB if (sum(abs(ipjpk(:)-int(ipjpk(:)))).lt.tol4) then
845 : !FB write(Invar%stdout,'(a,i4,1x,i4,1x,i4)')'OK for qi+qj+qk=G',iq_bz,jq_bz,kq_bz
846 : !FB okp_count=okp_count+1
847 : !FB else if (sum(abs(imjpk(:)-int(imjpk(:)))).lt.tol4) then
848 : !FB write(Invar%stdout,'(a,i4,1x,i4,1x,i4)')'OK for qi-qj+qk=G',iq_bz,jq_bz,kq_bz
849 : !FB okm_count=okm_count+1
850 : !FB else
851 : !FB cycle
852 : !FB end if
853 : !FB! do imode=1,nmode
854 : !FB! do jmode=1,nmode
855 : !FB! do kmode=1,nmode
856 : !FB! end do
857 : !FB! end do
858 : !FB! end do
859 : !FB end do
860 : !FB end do
861 : !FB end do
862 : !FB write(Invar%stdout,'(a,i10)') 'Total number of (q1,q2,q3) =',tot_count
863 : !FB write(Invar%stdout,'(a,i10)') 'Partial number of (q1+q2+q3) =',okp_count
864 : !FB write(Invar%stdout,'(a,i10)') 'Partial number of (q1-q2+q3) =',okm_count
865 : !FB
866 : !FBend subroutine tdep_calc_lifetime1
867 :
868 : !=====================================================================================================
869 : !FBsubroutine tdep_calc_lifetime2(Crystal,distance,Eigen2nd,Ifc,Invar,Lattice,MPIdata,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym)
870 : !FB
871 : !FB!Arguments ------------------------------------
872 : !FB!scalars
873 : !FB integer,intent(in) :: ncid,comm
874 : !FB real(dp),intent(in) :: dosdeltae !,dossmear
875 : !FB type(gruns_t),intent(in) :: gruns
876 : !FB character(len=*),intent(in) :: prefix
877 : !FB!arrays
878 : !FB
879 : !FB type(crystal_t),intent(in) :: Crystal
880 : !FB type(Symmetries_type),intent(in) :: Sym
881 : !FB type(atdep_dataset_type),intent(in) :: Invar
882 : !FB type(Shell_type),intent(in) :: Shell3at
883 : !FB type(Lattice_type),intent(in) :: Lattice
884 : !FB type(Eigen_type),intent(in) :: Eigen2nd
885 : !FB type(MPI_enreg_type), intent(in) :: MPIdata
886 : !FB type(Qbz_type),intent(in) :: Qbz
887 : !FB type(ifc_type),intent(in) :: Ifc
888 : !FB double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
889 : !FB double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
890 : !FB double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
891 : !FB
892 : !FB!Local variables-------------------------------
893 : !FB!scalars
894 : !FB integer,parameter :: master=0,qptopt1=1,bcorr0=0
895 : !FB integer :: nprocs,my_rank,iqibz,nqbz,nqibz,ierr,ii,nu,ncerr,nomega,cnt,unt,io
896 : !FB real(dp) :: gavg,omega_min,omega_max,v2
897 : !FB type(htetra_t) :: tetra
898 : !FB character(len=500) :: msg
899 : !FB!arrays
900 : !FB integer :: qptrlatt(3,3)
901 : !FB integer :: ngqpt(3)
902 : !FB real(dp) :: shiftq(3,1)
903 : !FB real(dp),allocatable :: gvals_qibz(:,:),wvols_qibz(:,:,:),dwdq_qibz(:,:,:)
904 : !FB real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:)
905 : !FB real(dp),allocatable :: wdt(:,:),wdos(:,:),grdos(:,:),gr2dos(:,:),wibz(:),omega_mesh(:)
906 : !FB real(dp),allocatable :: vdos(:,:),v2dos(:,:)
907 : !FB real(dp),allocatable :: phdispl_cart_qibz(:,:,:,:,:)
908 : !FB
909 : !FB! ************************************************************************
910 : !FB
911 : !FB comm = MPIdata%comm_step
912 : !FB nprocs = xmpi_comm_size(comm)
913 : !FB my_rank = xmpi_comm_rank(comm)
914 : !FB
915 : !FB write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Lifetime ',ch10
916 : !FB call wrtout(Invar%stdout,msg)
917 : !FB
918 : !FB ! Generate the q-mesh by finding the IBZ and the corresponding weights.
919 : !FB qptrlatt = 0
920 : !FB do ii=1,3
921 : !FB qptrlatt(ii,ii) = Invar%ngqpt1(:)
922 : !FB end do
923 : !FB nshiftq=1
924 : !FB shiftq(:,:)=0.0
925 : !FB
926 : !FB ! Get IBZ and BZ.
927 : !FB call kpts_ibz_from_kptrlatt(Crystal, qptrlatt, qptopt1, nshiftq, shiftq, &
928 : !FB nqibz, qibz, wtq, nqbz, qbz)
929 : !FB
930 : !FB ! Build tetrahedra
931 : !FB tetra = tetra_from_kptrlatt(Crystal, qptopt1, qptrlatt, nshiftq, shiftq, nqibz, qibz, comm, msg, ierr)
932 : !FB if (ierr /= 0) ABI_ERROR(msg)
933 : !FB
934 : !FB ABI_CALLOC(wvols_qibz, (gruns%natom3, gruns%nvols, nqibz))
935 : !FB ABI_CALLOC(gvals_qibz, (gruns%natom3, nqibz))
936 : !FB ABI_CALLOC(dwdq_qibz, (3, gruns%natom3, nqibz))
937 : !FB ABI_CALLOC(phdispl_cart_qibz, (2, gruns%natom3, gruns%natom3, gruns%nvols, nqibz))
938 : !FB
939 : !FB gavg = zero
940 : !FB do iqibz=1,nqibz
941 : !FB if (mod(iqibz, nprocs) /= my_rank) cycle ! mpi-parallelism
942 : !FB call gruns_fourq(gruns, qibz(:,iqibz), wvols_qibz(:,:,iqibz), gvals_qibz(:,iqibz), &
943 : !FB dwdq_qibz(:,:,iqibz), phdispl_cart_qibz(:,:,:,:,iqibz))
944 : !FB gavg = gavg + wtq(iqibz) * sum(gvals_qibz(:,iqibz))
945 : !FB end do
946 : !FB gavg = gavg / gruns%natom3
947 : !FB
948 : !FB call xmpi_sum(gavg, comm, ierr)
949 : !FB call xmpi_sum(wvols_qibz, comm, ierr)
950 : !FB call xmpi_sum(gvals_qibz, comm, ierr)
951 : !FB call xmpi_sum(dwdq_qibz, comm, ierr)
952 : !FB call xmpi_sum(phdispl_cart_qibz, comm, ierr)
953 : !FB
954 : !FB omega_min = gruns%ifc_vol(gruns%iv0)%omega_minmax(1)
955 : !FB omega_max = gruns%ifc_vol(gruns%iv0)%omega_minmax(2)
956 : !FB nomega = nint((omega_max - omega_min) / dosdeltae) + 1
957 : !FB nomega = max(6, nomega) ! Ensure Simpson integration will be ok
958 : !FB
959 : !FB ABI_MALLOC(omega_mesh, (nomega))
960 : !FB omega_mesh = arth(omega_min, dosdeltae, nomega)
961 : !FB omega_max = omega_mesh(nomega)
962 : !FB !write(std_out,*)"hello",omega_min,omega_max,dosdeltae,(omega_max-omega_min) / (nomega-1)
963 : !FB ABI_MALLOC(wibz, (nqibz))
964 : !FB ABI_MALLOC(wdt, (nomega, 2))
965 : !FB ABI_CALLOC(wdos, (nomega, 2))
966 : !FB ABI_CALLOC(grdos, (nomega, 2))
967 : !FB ABI_CALLOC(gr2dos, (nomega, 2))
968 : !FB ABI_CALLOC(vdos, (nomega, 2))
969 : !FB ABI_CALLOC(v2dos, (nomega, 2))
970 : !FB
971 : !FB ! Compute DOSes.
972 : !FB cnt = 0
973 : !FB do iqibz=1,nqibz
974 : !FB do nu=1,gruns%natom3
975 : !FB cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi-parallelism
976 : !FB wibz = wvols_qibz(nu, gruns%iv0, :)
977 : !FB call tetra%get_onewk(iqibz,bcorr0,nomega,nqibz,wibz,omega_min,omega_max,one,wdt)
978 : !FB wdt = wdt*wtq(iqibz)
979 : !FB wdos = wdos + wdt
980 : !FB grdos = grdos + wdt * gvals_qibz(nu,iqibz)
981 : !FB gr2dos = gr2dos + wdt * gvals_qibz(nu,iqibz) ** 2
982 : !FB v2 = sum(dwdq_qibz(1:3,nu,iqibz) ** 2)
983 : !FB vdos = vdos + wdt * sqrt(v2)
984 : !FB v2dos = v2dos + wdt * v2
985 : !FB end do
986 : !FB end do
987 : !FB
988 : !FB call xmpi_sum(wdos, comm, ierr)
989 : !FB call xmpi_sum(grdos, comm, ierr)
990 : !FB call xmpi_sum(gr2dos, comm, ierr)
991 : !FB call xmpi_sum(vdos, comm, ierr)
992 : !FB call xmpi_sum(v2dos, comm, ierr)
993 : !FB
994 : !FB if (my_rank == master) then
995 : !FB write(unt,'(a)')'# Phonon density of states, Gruneisen DOS and phonon group velocity DOS'
996 : !FB write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
997 : !FB write(unt,'(a,i0)')'# Tetrahedron method with nqibz= ',nqibz
998 : !FB write(unt,"(a,f8.5)")"# Average Gruneisen parameter:", gavg
999 : !FB write(unt,'(5a)') &
1000 : !FB "# omega PH_DOS Gruns_DOS Gruns**2_DOS Vel_DOS Vel**2_DOS PH_IDOS Gruns_IDOS Gruns**2_IDOS Vel_IDOS Vel**2_IDOS"
1001 : !FB do io=1,nomega
1002 : !FB write(unt, "(11es17.8)")omega_mesh(io), &
1003 : !FB wdos(io,1), grdos(io,1), gr2dos(io,1), vdos(io,1), v2dos(io,1), &
1004 : !FB wdos(io,2), grdos(io,2), gr2dos(io,2), vdos(io,2), v2dos(io,2)
1005 : !FB end do
1006 : !FB end if
1007 : !FB
1008 : !FB
1009 : !FB ABI_FREE(qibz)
1010 : !FB ABI_FREE(wtq)
1011 : !FB ABI_FREE(qbz)
1012 : !FB ABI_FREE(wvols_qibz)
1013 : !FB ABI_FREE(gvals_qibz)
1014 : !FB ABI_FREE(dwdq_qibz)
1015 : !FB ABI_FREE(phdispl_cart_qibz)
1016 : !FB ABI_FREE(omega_mesh)
1017 : !FB ABI_FREE(wibz)
1018 : !FB ABI_FREE(wdt)
1019 : !FB ABI_FREE(wdos)
1020 : !FB ABI_FREE(grdos)
1021 : !FB ABI_FREE(gr2dos)
1022 : !FB ABI_FREE(v2dos)
1023 : !FB ABI_FREE(vdos)
1024 : !FB
1025 : !FB call tetra%free()
1026 : !FB
1027 : !FBend subroutine tdep_calc_lifetime2
1028 : !!***
|
