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