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_phdos
9 :
10 : use defs_basis
11 : !FB use m_nctk
12 : use m_errors
13 : use m_abicore
14 : use m_phonons
15 : use m_xmpi
16 : use m_io_tools
17 : use m_ifc, only : ifc_type
18 : use m_crystal, only : crystal_t
19 : use m_ddb, only : ddb_type
20 : use m_tdep_phi2, only : Eigen_type, tdep_write_yaml, tdep_write_dij, tdep_calc_dij
21 : use m_tdep_model, only : tdep_init_phi2, tdep_destroy_phi2, Phi2_type
22 : use m_tdep_qpt, only : Qpoints_type
23 : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
24 : use m_tdep_latt, only : Lattice_type
25 : use m_tdep_sym, only : Symmetries_type
26 : use m_tdep_shell, only : Shell_type
27 : use m_tdep_abitypes, only : Qbz_type, tdep_ifc2phi2, tdep_read_ifc, tdep_write_ifc, &
28 : & tdep_write_ddb,tdep_init_ifc
29 :
30 : implicit none
31 :
32 : public :: tdep_calc_phdos
33 : public :: tdep_calc_thermo
34 : public :: tdep_calc_elastic
35 :
36 : contains
37 :
38 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
39 616 : subroutine tdep_calc_phdos(Crystal,DDB,Eigen2nd_MP,Eigen2nd_path,Ifc,Invar,Lattice,MPIdata,natom,&
40 44 : & natom_unitcell,Phi2,PHdos,Qbz,Qpt,Rlatt_scaled,Shell2at,Sym)
41 :
42 : integer, intent(in) :: natom,natom_unitcell
43 : double precision, intent(in) :: Rlatt_scaled(3,natom_unitcell,natom)
44 : type(atdep_dataset_type),intent(in) :: Invar
45 : type(phdos_t),intent(out) :: PHdos
46 : type(Phi2_type),intent(in) :: Phi2
47 : type(ifc_type),intent(inout) :: Ifc
48 : type(Lattice_type),intent(in) :: Lattice
49 : type(Symmetries_type),intent(in) :: Sym
50 : type(crystal_t),intent(inout) :: Crystal
51 : type(ddb_type),intent(in) :: DDB ! To be removed
52 44 : type(Phi2_type) :: Phi2_tmp ! To be removed
53 : type(Qbz_type),intent(in) :: Qbz
54 : type(Qpoints_type),intent(in) :: Qpt
55 : type(MPI_enreg_type), intent(in) :: MPIdata
56 : type(Shell_type),intent(in) :: Shell2at
57 : type(Eigen_type),intent(inout) :: Eigen2nd_path
58 : type(Eigen_type),intent(inout) :: Eigen2nd_MP
59 :
60 : integer :: prtdos,iqpt,iq_ibz,iomega,iatom
61 : integer :: dos_ngqpt(3)
62 : integer :: count_wminmax(2)
63 : character (len=fnlen):: phdos_fname
64 : double precision :: dossmear,integ,domega
65 : double precision :: dos_qshift(3)
66 44 : double precision, allocatable :: displ(:,:)
67 : character(len=500) :: message
68 : real(dp) :: wminmax(2)
69 264 : type(ifc_type) :: Ifc_tmp
70 :
71 : !FB double precision, allocatable :: omega (:)
72 : !FB double complex, allocatable :: eigenV(:,:)
73 :
74 44 : write(Invar%stdout,*)' '
75 44 : write(Invar%stdout,*) '#############################################################################'
76 44 : write(Invar%stdout,*) '############## Compute the phonon spectrum, the DOS, ########################'
77 44 : write(Invar%stdout,*) '############## the dynamical matrix and write them ########################'
78 44 : write(Invar%stdout,*) '#############################################################################'
79 :
80 44 : write(Invar%stdout,'(a)')' '
81 44 : write(Invar%stdout,'(a)') ' #############################################################################'
82 44 : write(Invar%stdout,'(a)') ' ################### vibrational Density OF States (vDOS) ####################'
83 44 : write(Invar%stdout,'(a)') ' #############################################################################'
84 44 : write(Invar%stdout,'(a)') ' See the vdos.dat and TDEP_PHDOS* files'
85 :
86 : ! ----------------------------------------------------
87 : ! GA To be moved outside of this routine vvvvvvvvvvvvv
88 :
89 : ! Write Ifc%atmfrc in the ifc_out.dat file
90 44 : if (MPIdata%iam_master) then
91 44 : call tdep_write_ifc(Crystal,Ifc,Invar,Invar%natom_unitcell,0)
92 : end if
93 :
94 : ! For test purpose : read the previous IFC from ifc_out.dat and write it in ifc_check.dat
95 : ! =======================================================================================
96 44 : if (Invar%readifc.eq.2) then
97 :
98 : ! Copy Phi2 to Phi2_tmp
99 2 : call tdep_init_phi2(Phi2_tmp,Invar%loto,natom)
100 253144 : Phi2_tmp%SR =Phi2%SR
101 2 : if (Invar%loto) then
102 0 : Phi2_tmp%Tot=Phi2%Tot
103 0 : Phi2_tmp%LR =Phi2%LR
104 : end if
105 :
106 : ! Copy Phi2_tmp to Ifc_tmp%atmfrc
107 2 : call tdep_init_ifc(Crystal,DDB,Ifc_tmp,Invar,Lattice,MPIdata,Phi2_tmp,Rlatt_scaled,Shell2at,Sym)
108 :
109 2 : if (MPIdata%iam_master) then
110 :
111 : ! Read IFC from ifc_out.dat (readifc=2)
112 2 : call tdep_read_ifc(Ifc_tmp,Invar,natom_unitcell)
113 :
114 : ! Copy Ifc_tmp%atmfrc to Phi2_tmp
115 2 : call tdep_ifc2phi2(Ifc_tmp%dipdip,Ifc_tmp,Invar,Lattice,natom_unitcell,1,Phi2_tmp,Rlatt_scaled,Shell2at,Sym)
116 :
117 : ! Copy Phi2_tmp to Ifc_tmp%atmfrc
118 2 : call tdep_ifc2phi2(Ifc_tmp%dipdip,Ifc_tmp,Invar,Lattice,natom_unitcell,0,Phi2_tmp,Rlatt_scaled,Shell2at,Sym)
119 :
120 : ! Write IFC in ifc_check.dat (for check)
121 2 : call tdep_write_ifc(Crystal,Ifc_tmp,Invar,natom_unitcell,1)
122 :
123 : ! Write the Phi2-tmp.dat file
124 2 : if (Invar%debug) then
125 0 : write(Invar%stdout,'(a)') ' See the Phi2-tmp.dat file corresponding to the ifc_out.dat/Phi2 file'
126 0 : open(unit=55,file=trim(Invar%output_prefix)//'_Phi2-tmp.dat')
127 0 : do iatom=1,3*natom
128 0 : write(55,'(10000(f10.6,1x))') Phi2_tmp%SR(iatom,:)
129 : end do
130 0 : close(55)
131 : end if
132 : end if
133 2 : call Ifc_tmp%free()
134 2 : call tdep_destroy_phi2(Phi2_tmp)
135 : end if
136 :
137 : ! GA To be moved outside of this routine ^^^^^^^^^^^^^
138 : ! ----------------------------------------------------
139 :
140 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ON THE FINE GRID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143 : ! Compute the DOS
144 : ! ===============
145 44 : prtdos=1 !Gaussian
146 : ! prtdos=2 !Tetra
147 44 : dossmear=4.5d-6
148 : ! dossmear=4.5d-5
149 176 : dos_qshift(:)=0.5d0
150 176 : dos_ngqpt(:)=Invar%ngqpt2(:)
151 44 : write(Invar%stdout,'(a)') ' Compute the vDOS'
152 : ! Only 1 shift in q-mesh
153 44 : wminmax = zero
154 44 : do
155 : call PHdos%init(Crystal,Ifc,prtdos,Invar%dosdeltae,dossmear,dos_ngqpt,1,dos_qshift, &
156 88 : "freq_displ", wminmax, count_wminmax, XMPI_WORLD,prtout=.false.)
157 220 : if (all(count_wminmax == 0)) exit
158 44 : wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05
159 44 : wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05
160 44 : call phdos%free()
161 44 : write(message, "(a, 2f8.5)")"Initial frequency mesh not large enough. Recomputing PHDOS with wmin, wmax: ",wminmax
162 44 : call wrtout(std_out, message)
163 : end do
164 44 : write(Invar%stdout,'(a)') ' ------- achieved'
165 44 : write(Invar%stdout,'(a)') ' (Please, pay attention to convergency wrt the BZ mesh : the ngqpt2 input variable)'
166 :
167 : ! Print the DOS
168 : ! =============
169 44 : phdos_fname = trim(Invar%output_prefix)//"_PHDOS"
170 44 : if (MPIdata%iam_master) call phdos%print(phdos_fname)
171 44 : domega=(Invar%dosdeltae*Ha_meV)
172 44 : integ=0.d0
173 80267 : do iomega=1,PHdos%nomega
174 80267 : integ=integ + domega*PHdos%phdos(iomega)
175 : end do
176 80267 : PHdos%phdos(:)=PHdos%phdos(:)/integ
177 44 : if (MPIdata%iam_master) then
178 44 : open(unit=56,file=trim(Invar%output_prefix)//'_vdos.dat')
179 80267 : do iomega=1,PHdos%nomega
180 80223 : if (Invar%enunit.eq.0) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_eV*1000,PHdos%phdos(iomega)
181 80223 : if (Invar%enunit.eq.1) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_cmm1 ,PHdos%phdos(iomega)
182 80223 : if (Invar%enunit.eq.2) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*1000 ,PHdos%phdos(iomega)
183 80267 : if (Invar%enunit.eq.3) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_THz ,PHdos%phdos(iomega)
184 : end do
185 44 : close(56)
186 : end if
187 :
188 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ON THE PATH GRID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191 : ! Compute the spectrum and the dynamical matrix
192 : ! =============================================
193 : ! Compute the frequencies (1)
194 : ! =======================
195 : !FB Eigen2nd_path%eigenval=zero ; Eigen2nd_path%eigenvec=zero ; Eigen2nd_path%dynmat=zero
196 : !FB ABI_MALLOC(dij ,(3*Invar%natom_unitcell,3*Invar%natom_unitcell))
197 : !FB ABI_MALLOC(eigenV,(3*Invar%natom_unitcell,3*Invar%natom_unitcell))
198 : !FB ABI_MALLOC(omega,(3*Invar%natom_unitcell))
199 : !FB do iqpt=1,Qpt%nqpt
200 : !FB omega=zero ; eigenV=zero ; dij=zero
201 : !FB call tdep_calc_dij(dij,eigenV,iqpt,Invar,omega,Phi2%Tot,Qpt%qpt_cart(:,iqpt),Rlatt_cart)
202 : !FB Eigen2nd_path%eigenval(:,iqpt)= omega(:)
203 : !FB do iatom=1,Invar%natom_unitcell
204 : !FB do ii=1,3
205 : !FB do jatom=1,Invar%natom_unitcell
206 : !FB do jj=1,3
207 : !FB Eigen2nd_path%eigenvec(1,ii,iatom,jj,jatom,iqpt)= real(eigenV(ii+3*(iatom-1),jj+3*(jatom-1)))
208 : !FB Eigen2nd_path%eigenvec(2,ii,iatom,jj,jatom,iqpt)=aimag(eigenV(ii+3*(iatom-1),jj+3*(jatom-1)))
209 : !FB Eigen2nd_path%dynmat(1,ii,iatom,jj,jatom,iqpt) = real(dij(ii+3*(iatom-1),jj+3*(jatom-1)))
210 : !FB Eigen2nd_path%dynmat(2,ii,iatom,jj,jatom,iqpt) =aimag(dij(ii+3*(iatom-1),jj+3*(jatom-1)))
211 : !FB end do
212 : !FB end do
213 : !FB end do
214 : !FB end do
215 : !FB end do
216 : !FB ABI_FREE(dij)
217 : !FB ABI_FREE(eigenV)
218 : !FB ABI_FREE(omega)
219 : !FB! Write the Dij, eigenvalues and eigenvectors (in ASCII)
220 : !FB! ======================================================
221 : !FB if (MPIdata%iam_master) then
222 : !FB open(unit=51,file=trim(Invar%output_prefix)//'eigenvectors-path-1.dat')
223 : !FB open(unit=52,file=trim(Invar%output_prefix)//'dij-path-1.dat')
224 : !FB open(unit=53,file=trim(Invar%output_prefix)//'omega-path-1.dat')
225 : !FB if (Invar%enunit.eq.0) write(53,'(a)') '# Phonon frequencies in meV'
226 : !FB if (Invar%enunit.eq.1) write(53,'(a)') '# Phonon frequencies in cm-1'
227 : !FB if (Invar%enunit.eq.2) write(53,'(a)') '# Phonon frequencies in mHa'
228 : !FB if (Invar%enunit.eq.3) write(53,'(a)') '# Phonon frequencies in THz'
229 : !FB do iqpt=1,Qpt%nqpt
230 : !FB call tdep_write_dij(Eigen2nd_path,iqpt,Invar,Qpt%qpt_red(:,iqpt))
231 : !FB end do
232 : !FB close(51)
233 : !FB close(52)
234 : !FB close(53)
235 : !FB end if
236 :
237 : ! Compute the frequencies (2)
238 : ! =======================
239 3642042 : ABI_MALLOC(displ,(2*3*natom_unitcell*3*natom_unitcell,Qpt%nqpt)); displ(:,:)=zero
240 12694465 : Eigen2nd_path%eigenval=zero ; Eigen2nd_path%eigenvec=zero ; Eigen2nd_path%dynmat=zero
241 22632 : do iqpt=1,Qpt%nqpt
242 : call Ifc%fourq(Crystal,Qpt%qpt_red(:,iqpt),Eigen2nd_path%eigenval(:,iqpt),displ(:,iqpt),&
243 : & out_eigvec=Eigen2nd_path%eigenvec(:,:,:,:,:,iqpt),&
244 22632 : & out_d2cart=Eigen2nd_path%dynmat (:,:,:,:,:,iqpt))
245 : end do
246 44 : ABI_FREE(displ)
247 : ! Write the Dij, eigenvalues and eigenvectors (in ASCII and YAML)
248 : ! ======================================================
249 44 : if (MPIdata%iam_master) then
250 : !FB open(unit=51,file=trim(Invar%output_prefix)//'eigenvectors-path-2.dat')
251 : !FB open(unit=52,file=trim(Invar%output_prefix)//'dij-path-2.dat')
252 : !FB open(unit=53,file=trim(Invar%output_prefix)//'omega-path-2.dat')
253 44 : open(unit=51,file=trim(Invar%output_prefix)//'_eigenvectors.dat')
254 44 : open(unit=52,file=trim(Invar%output_prefix)//'_dij.dat')
255 44 : open(unit=53,file=trim(Invar%output_prefix)//'_omega.dat')
256 44 : if (Invar%enunit.eq.0) write(53,'(a)') '# Phonon frequencies in meV'
257 44 : if (Invar%enunit.eq.1) write(53,'(a)') '# Phonon frequencies in cm-1'
258 44 : if (Invar%enunit.eq.2) write(53,'(a)') '# Phonon frequencies in mHa'
259 44 : if (Invar%enunit.eq.3) write(53,'(a)') '# Phonon frequencies in THz'
260 22632 : do iqpt=1,Qpt%nqpt
261 22632 : call tdep_write_dij(Eigen2nd_path,iqpt,Invar,Qpt%qpt_red(:,iqpt))
262 : end do
263 44 : close(51)
264 44 : close(52)
265 44 : close(53)
266 44 : call tdep_write_yaml(Eigen2nd_path,Qpt,Invar%output_prefix)
267 : end if
268 :
269 : !FB!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
270 : !FB!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ON THE MP GRID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
271 : !FB!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272 : !FB! Compute the spectrum and the dynamical matrix
273 : !FB! =============================================
274 : !FB if (Qbz%nqbz.ne.Qbz%nqibz) then
275 : !FB write(Invar%stdlog,*) 'DIFFERENCES BETWEEN NQBZ \& NQIBZ :',Qbz%nqbz,Qbz%nqibz
276 : !FB end if
277 : !FB! Compute the frequencies(1)
278 : !FB! =======================
279 : !FB Eigen2nd_MP%eigenval=zero ; Eigen2nd_MP%eigenvec=zero ; Eigen2nd_MP%dynmat=zero
280 : !FB ABI_MALLOC(dij ,(3*Invar%natom_unitcell,3*Invar%natom_unitcell))
281 : !FB ABI_MALLOC(eigenV,(3*Invar%natom_unitcell,3*Invar%natom_unitcell))
282 : !FB ABI_MALLOC(omega,(3*Invar%natom_unitcell))
283 : !FB do iq_ibz=1,Qbz%nqibz
284 : !FB omega=zero ; eigenV=zero ; dij=zero
285 : !FB call tdep_calc_dij(dij,eigenV,iq_ibz,Invar,omega,Phi2,Qbz%qibz_cart(:,iq_ibz),Rlatt_cart)
286 : !FB Eigen2nd_MP%eigenval(:,iq_ibz)= omega(:)
287 : !FB do iatom=1,Invar%natom_unitcell
288 : !FB do ii=1,3
289 : !FB do jatom=1,Invar%natom_unitcell
290 : !FB do jj=1,3
291 : !FB Eigen2nd_MP%eigenvec(1,ii,iatom,jj,jatom,iq_ibz)= real(eigenV(ii+3*(iatom-1),jj+3*(jatom-1)))
292 : !FB Eigen2nd_MP%eigenvec(2,ii,iatom,jj,jatom,iq_ibz)=aimag(eigenV(ii+3*(iatom-1),jj+3*(jatom-1)))
293 : !FB Eigen2nd_MP%dynmat(1,ii,iatom,jj,jatom,iq_ibz) = real(dij(ii+3*(iatom-1),jj+3*(jatom-1)))
294 : !FB Eigen2nd_MP%dynmat(2,ii,iatom,jj,jatom,iq_ibz) =aimag(dij(ii+3*(iatom-1),jj+3*(jatom-1)))
295 : !FB end do
296 : !FB end do
297 : !FB end do
298 : !FB end do
299 : !FB end do
300 : !FB ABI_FREE(dij)
301 : !FB ABI_FREE(eigenV)
302 : !FB ABI_FREE(omega)
303 : !FB if (MPIdata%iam_master) then
304 : !FB open(unit=51,file=trim(Invar%output_prefix)//'eigenvectors-MP-1.dat')
305 : !FB open(unit=52,file=trim(Invar%output_prefix)//'dij-MP-1.dat')
306 : !FB open(unit=53,file=trim(Invar%output_prefix)//'omega-MP-1.dat')
307 : !FB if (Invar%enunit.eq.0) write(53,'(a)') '# Phonon frequencies in meV'
308 : !FB if (Invar%enunit.eq.1) write(53,'(a)') '# Phonon frequencies in cm-1'
309 : !FB if (Invar%enunit.eq.2) write(53,'(a)') '# Phonon frequencies in mHa'
310 : !FB if (Invar%enunit.eq.3) write(53,'(a)') '# Phonon frequencies in THz'
311 : !FB do iq_ibz=1,Qbz%nqibz
312 : !FB call tdep_write_dij(Eigen2nd_MP,iq_ibz,Invar,Qbz%qibz(:,iq_ibz))
313 : !FB end do
314 : !FB close(51)
315 : !FB close(52)
316 : !FB close(53)
317 : !FB end if
318 : !FB
319 : ! Compute the frequencies (2)
320 : ! =======================
321 2109608 : ABI_MALLOC(displ,(2*3*natom_unitcell*3*natom_unitcell,Qbz%nqibz)); displ(:,:)=zero
322 14478780 : Eigen2nd_MP%eigenval=zero ; Eigen2nd_MP%eigenvec=zero ; Eigen2nd_MP%dynmat=zero
323 7076 : do iq_ibz=1,Qbz%nqibz
324 : call Ifc%fourq(Crystal,Qbz%qibz(:,iq_ibz),Eigen2nd_MP%eigenval(:,iq_ibz),displ(:,iq_ibz),&
325 : & out_eigvec=Eigen2nd_MP%eigenvec(:,:,:,:,:,iq_ibz),&
326 7076 : & out_d2cart=Eigen2nd_MP%dynmat (:,:,:,:,:,iq_ibz))
327 : end do
328 44 : ABI_FREE(displ)
329 : !FB if (MPIdata%iam_master) then
330 : !FB open(unit=51,file=trim(Invar%output_prefix)//'eigenvectors-MP-2.dat')
331 : !FB open(unit=52,file=trim(Invar%output_prefix)//'dij-MP-2.dat')
332 : !FB open(unit=53,file=trim(Invar%output_prefix)//'omega-MP-2.dat')
333 : !FB if (Invar%enunit.eq.0) write(53,'(a)') '# Phonon frequencies in meV'
334 : !FB if (Invar%enunit.eq.1) write(53,'(a)') '# Phonon frequencies in cm-1'
335 : !FB if (Invar%enunit.eq.2) write(53,'(a)') '# Phonon frequencies in mHa'
336 : !FB if (Invar%enunit.eq.3) write(53,'(a)') '# Phonon frequencies in THz'
337 : !FB do iq_ibz=1,Qbz%nqibz
338 : !FB call tdep_write_dij(Eigen2nd_MP,iq_ibz,Invar,Qbz%qibz(:,iq_ibz))
339 : !FB end do
340 : !FB close(51)
341 : !FB close(52)
342 : !FB close(53)
343 : !FB end if
344 :
345 44 : end subroutine tdep_calc_phdos
346 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347 :
348 44 : subroutine tdep_calc_thermo(Invar,Lattice,MPIdata,PHdos,U0)
349 :
350 : double precision, intent(in) :: U0
351 : type(atdep_dataset_type),intent(in) :: Invar
352 : type(Lattice_type), intent(inout) :: Lattice
353 : type(MPI_enreg_type), intent(in) :: MPIdata
354 : type(phdos_t),intent(in) :: PHdos
355 :
356 : integer :: iomega,itemp,iatom,itypat
357 : double precision :: k_B,wovert,heatcapa,entropy,internalE,freeE,expm2x,ln2shx,cothx,xx
358 : double precision :: Ftot,domega,MSD,Omega_m2,mass_amu,vdos
359 :
360 44 : write(Invar%stdout,'(a)')' '
361 44 : write(Invar%stdout,'(a)') ' #############################################################################'
362 44 : write(Invar%stdout,'(a)') ' ################# Thermodynamic quantities: Free energy,...##################'
363 44 : write(Invar%stdout,'(a)') ' #############################################################################'
364 44 : write(Invar%stdout,'(a)') ' See the thermo.dat file'
365 :
366 : ! The heat capacity, entropy, internal and free energies (direct calculation)
367 : ! ===========================================================================
368 44 : mass_amu=zero
369 150 : do iatom=1,Invar%natom_unitcell
370 106 : itypat=Invar%typat_unitcell(iatom)
371 150 : mass_amu=mass_amu+Invar%amu(itypat)
372 : end do
373 44 : mass_amu=mass_amu*amu_emass/real(Invar%natom_unitcell)
374 :
375 : ! k_B is Boltzmann constant in eV/K
376 44 : k_B=kb_HaK*Ha_eV
377 44 : domega=(Invar%dosdeltae*Ha_meV)
378 44 : wovert=1.d0/(2*Invar%temperature*k_B)
379 44 : heatcapa=0.d0 ; entropy=0.d0 ; internalE=0.d0 ; freeE=0.d0 ; MSD=0.d0 ; Omega_m2=0.d0 ; vdos=0.d0
380 80267 : do iomega=1,PHdos%nomega
381 80223 : xx=PHdos%omega(iomega)*Ha_eV
382 80223 : if (xx.lt.tol8) cycle
383 80015 : expm2x=exp(-2.d0*wovert*xx)
384 80015 : ln2shx=wovert*xx+log(1.d0-expm2x)
385 80015 : cothx=(1.d0+expm2x)/(1.d0-expm2x)
386 80015 : heatcapa =heatcapa + (wovert*xx/sinh(wovert*xx))**2*PHdos%phdos(iomega)*domega
387 80015 : internalE=internalE + (wovert*xx)*cothx*PHdos%phdos(iomega)*domega
388 80015 : entropy =entropy + ((wovert*xx)*cothx-ln2shx)*PHdos%phdos(iomega)*domega
389 80015 : freeE =freeE + log(2*sinh(wovert*xx))*PHdos%phdos(iomega)*domega
390 80015 : MSD =MSD + (cothx/PHdos%omega(iomega))*PHdos%phdos(iomega)*domega
391 80015 : Omega_m2 =Omega_m2 + (1.d0/PHdos%omega(iomega))**2*PHdos%phdos(iomega)*domega
392 80267 : vdos =vdos + PHdos%phdos(iomega)*domega
393 : end do
394 44 : heatcapa=heatcapa*3
395 44 : entropy=entropy*3
396 44 : internalE=internalE*3*k_B*Invar%temperature
397 44 : freeE=freeE*3*k_B*Invar%temperature
398 44 : MSD=MSD*3.d0/mass_amu/2.d0
399 44 : Omega_m2=Omega_m2*3.d0
400 44 : Ftot=U0*Ha_eV+freeE
401 44 : Lattice%HeatCapa_V=heatcapa
402 44 : if (.not.MPIdata%iam_master) return
403 : ! End of the calculation --> RETURN
404 : ! =================================
405 :
406 44 : open(unit=20,file=trim(Invar%output_prefix)//'_thermo.dat')
407 44 : write(20,'(a)')'============= Direct results (without any inter/extrapolation) =================='
408 44 : write(20,'(1x,a,f10.3)')'For present temperature (in Kelvin): T= ',Invar%temperature
409 44 : write(20,'(1x,a,f12.3)')' The cold contribution (in eV/atom): U_0 =',U0*Ha_eV
410 44 : write(20,'(1x,a,f10.3)')' The specific heat (in k_b/atom): C_v=',heatcapa
411 44 : write(20,'(1x,a,f10.3)')' The vibrational entropy (in k_b/atom): S_vib =',entropy
412 44 : write(20,'(1x,a,f10.3)')' The internal energy (in eV/atom): U_vib =',internalE
413 44 : write(20,'(1x,a,f10.3)')' The vibrational contribution (in eV/atom): F_vib = U_vib -T.S_vib =',freeE
414 44 : write(20,'(1x,a,f12.3)')' The harmonic free energy (in eV/atom) --> F_tot^HA = U_0 + F_vib =',Ftot
415 44 : write(20,'(1x,a)')' Useful quantities for melting :'
416 44 : write(20,'(1x,a,f10.3)')' The mean square displacement (in a.u.): sqrt(<u^2>) =',(MSD)**0.5
417 44 : write(20,'(1x,a,f10.3)')' The <Omega^(-2)> factor (in THz^(-2)) =',Omega_m2/(Ha_THz)**2
418 44 : write(20,'(1x,a,f10.3)')' The Wigner-Seitz radius (in a.u.) : d_at =',(6*Lattice%ucvol/pi/real(Invar%natom_unitcell))**(1./3.)
419 44 : write(20,'(1x,a,f10.3)')' The average mass / proton-electron mass ratio (in a.u.) =', mass_amu/amu_emass
420 44 : write(20,'(1x,a,f10.3)')' The Lindemann constant : sqrt(<u^2>)/d_at =',(MSD)**0.5/&
421 88 : & (6*Lattice%ucvol/pi/real(Invar%natom_unitcell))**(1./3.)
422 44 : write(20,'(1x,a,f10.3)')' The integral of vDOS =',vdos
423 44 : write(20,'(a)')' '
424 :
425 : ! The free energy (extrapolation)
426 : ! ===============================
427 44 : write(20,'(a)')'============= Harmonic Approximation (HA) =================='
428 44 : write(20,'(1x,a)')' Note that the following results come from an EXTRAPOLATION:'
429 44 : write(20,'(1x,a,i5,a)')' 1/ F_vib^HA(T) is computed for each T using vDOS(T=',int(Invar%temperature),')'
430 44 : write(20,'(1x,a)')' 2/ F_tot^HA(T) = F_vib^HA(T) + U_0'
431 : write(20,'(a)')' T F_vib^HA(T) F_tot^HA(T) C_v(T) '&
432 44 : & //' S_vib(T) U_vib(T) MSD(T)'
433 4444 : do itemp=1,100
434 4400 : wovert=1.d0/(2*real(itemp)*100*k_B)
435 4400 : freeE=0.d0
436 4400 : heatcapa=0.d0
437 4400 : entropy=0.d0
438 4400 : internalE=0.d0
439 4400 : MSD=0.d0
440 8026700 : do iomega=1,PHdos%nomega
441 8022300 : xx=PHdos%omega(iomega)*Ha_eV
442 8022300 : if (xx.lt.tol8) cycle
443 8001500 : expm2x=exp(-2.d0*wovert*xx)
444 8001500 : ln2shx=wovert*xx+log(1.d0-expm2x)
445 8001500 : cothx=(1.d0+expm2x)/(1.d0-expm2x)
446 8001500 : heatcapa =heatcapa + (wovert*xx/sinh(wovert*xx))**2*PHdos%phdos(iomega)*domega
447 8001500 : internalE=internalE + (wovert*xx)*cothx*PHdos%phdos(iomega)*domega
448 8001500 : entropy =entropy + ((wovert*xx)*cothx-ln2shx)*PHdos%phdos(iomega)*domega
449 8001500 : freeE =freeE + log(2*sinh(wovert*xx))*PHdos%phdos(iomega)*domega
450 8026700 : MSD =MSD + (cothx/PHdos%omega(iomega))*PHdos%phdos(iomega)*domega
451 : end do
452 4400 : heatcapa=heatcapa*3
453 4400 : entropy=entropy*3
454 4400 : internalE=internalE*3*k_B*itemp*100
455 4400 : freeE=freeE*3*k_B*itemp*100
456 4400 : MSD=MSD*3.d0/mass_amu/2.d0
457 4400 : Ftot=U0*Ha_eV+freeE
458 4444 : write(20,'(1x,i5,6(1x,f15.3))') itemp*100,freeE,Ftot,heatcapa,entropy,internalE,(MSD)**0.5
459 : end do
460 44 : close(20)
461 :
462 : end subroutine tdep_calc_thermo
463 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464 :
465 44 : subroutine tdep_calc_elastic(Phi2,distance,Invar,Lattice)
466 :
467 : integer :: iatom,ii,jj,kk,ll,iatcell,itypat
468 : integer :: INFO,LWORK
469 : double precision :: BH,BR,BV,GR,GV,GH,Eaverage,Nuaverage,Laverage,Vp,Vs,Vphi
470 : double precision :: rho,E1,E2,E3,Nu12,Nu13,Nu23,Nu21,Nu31,Nu32,G23,G13,G12
471 : double precision :: mass_amu,bohr,V_Debye,T_Debye,A_U,A_B,A_G
472 44 : integer, allocatable :: IPIV(:)
473 44 : double precision, allocatable :: eigenvalues(:)
474 44 : double precision, allocatable :: WORK(:)
475 44 : double precision, allocatable :: Sij(:,:),Cij(:,:),aijkl(:,:,:,:),cijkl(:,:,:,:)
476 : type(atdep_dataset_type), intent(in) :: Invar
477 : type(Lattice_type), intent(inout) :: Lattice
478 : double precision, intent(in) :: distance(Invar%natom,Invar%natom,4)
479 : double precision, intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
480 :
481 44 : write(Invar%stdout,'(a)')' '
482 44 : write(Invar%stdout,'(a)') ' #############################################################################'
483 44 : write(Invar%stdout,'(a)') ' ######################### Elastic constants #################################'
484 44 : write(Invar%stdout,'(a)') ' ################ Bulk and Shear modulus--Sound velocities ###################'
485 44 : write(Invar%stdout,'(a)') ' #############################################################################'
486 :
487 44 : bohr=Bohr_Ang*1.0e-10
488 : ! Define atomic mass average
489 44 : mass_amu=zero
490 150 : do iatom=1,Invar%natom_unitcell
491 106 : itypat=Invar%typat_unitcell(iatom)
492 150 : mass_amu=mass_amu+Invar%amu(itypat)
493 : end do
494 44 : mass_amu=mass_amu/real(Invar%natom_unitcell)
495 :
496 44 : rho=(mass_amu/1e3)*Invar%natom_unitcell/Lattice%ucvol/bohr**3/Avogadro
497 :
498 : !==========================================================================================
499 : !===================== Elastic constants ==================================================
500 : !==========================================================================================
501 : ! New calculation of elastic constants using the formula (12.28 and 12.29 of
502 : ! Wallace, Statistical physics of crystals and liquids, Worl Scientific)
503 5324 : ABI_MALLOC(aijkl,(3,3,3,3)); aijkl(:,:,:,:)=0.d0
504 5324 : ABI_MALLOC(cijkl,(3,3,3,3)); cijkl(:,:,:,:)=0.d0
505 176 : do ii=1,3
506 572 : do jj=1,3
507 1716 : do kk=1,3
508 5148 : do ll=1,3
509 562194 : do iatom=1,Invar%natom
510 1766772 : do iatcell=1,Invar%natom_unitcell
511 : aijkl(ii,jj,kk,ll)=aijkl(ii,jj,kk,ll)-Phi2(ii+(iatcell-1)*3,3*(iatom-1)+jj)&
512 1763208 : & *distance(iatcell,iatom,kk+1)*distance(iatcell,iatom,ll+1)/2.d0/Lattice%ucvol
513 : end do
514 : end do
515 : enddo
516 : enddo
517 : enddo
518 : enddo
519 :
520 176 : do ii=1,3
521 572 : do jj=1,3
522 1716 : do kk=1,3
523 5148 : do ll=1,3
524 4752 : cijkl(ii,jj,kk,ll)=aijkl(ii,kk,jj,ll)+aijkl(jj,kk,ii,ll)-aijkl(ii,jj,kk,ll)
525 : !FB cijkl(ii,jj,kk,ll)=aijkl(ii,kk,jj,ll)+aijkl(jj,kk,ii,ll)-aijkl(jj,ii,kk,ll)
526 : enddo
527 : enddo
528 : enddo
529 : enddo
530 44 : ABI_FREE(aijkl)
531 :
532 5324 : cijkl(:,:,:,:)=cijkl(:,:,:,:)*29421.033d0
533 :
534 1892 : ABI_MALLOC(Cij,(6,6)) ; Cij(:,:)=0.d0
535 44 : Cij(1,1)=cijkl(1,1,1,1) ; Cij(1,2)=cijkl(1,1,2,2) ; Cij(1,3)=cijkl(1,1,3,3)
536 44 : Cij(1,4)=cijkl(1,1,2,3) ; Cij(1,5)=cijkl(1,1,1,3) ; Cij(1,6)=cijkl(1,1,1,2)
537 44 : Cij(2,1)=cijkl(2,2,1,1) ; Cij(2,2)=cijkl(2,2,2,2) ; Cij(2,3)=cijkl(2,2,3,3)
538 44 : Cij(2,4)=cijkl(2,2,2,3) ; Cij(2,5)=cijkl(2,2,1,3) ; Cij(2,6)=cijkl(2,2,1,2)
539 44 : Cij(3,1)=cijkl(3,3,1,1) ; Cij(3,2)=cijkl(3,3,2,2) ; Cij(3,3)=cijkl(3,3,3,3)
540 44 : Cij(3,4)=cijkl(3,3,2,3) ; Cij(3,5)=cijkl(3,3,1,3) ; Cij(3,6)=cijkl(3,3,1,2)
541 44 : Cij(4,1)=cijkl(2,3,1,1) ; Cij(4,2)=cijkl(2,3,2,2) ; Cij(4,3)=cijkl(2,3,3,3)
542 44 : Cij(4,4)=cijkl(2,3,2,3) ; Cij(4,5)=cijkl(2,3,1,3) ; Cij(4,6)=cijkl(2,3,1,2)
543 44 : Cij(5,1)=cijkl(1,3,1,1) ; Cij(5,2)=cijkl(1,3,2,2) ; Cij(5,3)=cijkl(1,3,3,3)
544 44 : Cij(5,4)=cijkl(1,3,2,3) ; Cij(5,5)=cijkl(1,3,1,3) ; Cij(5,6)=cijkl(1,3,1,2)
545 44 : Cij(6,1)=cijkl(1,2,1,1) ; Cij(6,2)=cijkl(1,2,2,2) ; Cij(6,3)=cijkl(1,2,3,3)
546 44 : Cij(6,4)=cijkl(1,2,2,3) ; Cij(6,5)=cijkl(1,2,1,3) ; Cij(6,6)=cijkl(1,2,1,2)
547 : ! Remove the rounding errors before writing (for non regression testing purposes)
548 308 : do ii=1,6
549 1892 : do jj=1,6
550 1848 : if (abs(Cij(ii,jj)).lt.tol8) Cij(ii,jj)=zero
551 : end do
552 : end do
553 44 : write(Invar%stdout,'(a)') ' '
554 44 : write(Invar%stdout,'(a)') ' ========== Using the formulation proposed by Wallace (using the IFC) ========='
555 44 : write(Invar%stdout,'(a)') ' Cijkl [in GPa]='
556 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C11 C12 C13 C14 C15 C16 | ',Cij(1,1),Cij(1,2),Cij(1,3),Cij(1,4),Cij(1,5),Cij(1,6)
557 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C21 C22 C23 C24 C25 C26 | ',Cij(2,1),Cij(2,2),Cij(2,3),Cij(2,4),Cij(2,5),Cij(2,6)
558 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C31 C32 C33 C34 C35 C36 | ',Cij(3,1),Cij(3,2),Cij(3,3),Cij(3,4),Cij(3,5),Cij(3,6)
559 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C41 C42 C43 C44 C45 C46 | = ',Cij(4,1),Cij(4,2),Cij(4,3),Cij(4,4),Cij(4,5),Cij(4,6)
560 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C51 C52 C53 C54 C55 C56 | ',Cij(5,1),Cij(5,2),Cij(5,3),Cij(5,4),Cij(5,5),Cij(5,6)
561 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | C61 C62 C63 C64 C65 C66 | ',Cij(6,1),Cij(6,2),Cij(6,3),Cij(6,4),Cij(6,5),Cij(6,6)
562 :
563 : ! Compute the eigenvalues of the Cij matrix in order to find the Born-Huang
564 : ! stability criterion (see Wallace, Thermodynamics of crystal, p39)
565 44 : ABI_MALLOC(WORK,(1))
566 1892 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
567 308 : ABI_MALLOC(eigenvalues,(6)) ; eigenvalues(:)=0.d0
568 1892 : Sij(:,:)=Cij(:,:)
569 44 : LWORK=-1
570 44 : call DSYEV('N','U',6,Sij,6,eigenvalues,WORK,LWORK,INFO)
571 44 : LWORK=WORK(1)
572 44 : ABI_FREE(WORK)
573 132 : ABI_MALLOC(WORK,(LWORK))
574 44 : call DSYEV('N','U',6,Sij,6,eigenvalues,WORK,LWORK,INFO)
575 308 : do ii=1,6
576 308 : if (eigenvalues(ii).lt.0.d0) then
577 0 : write(Invar%stdout,'(a)') ' WARNING :'
578 0 : write(Invar%stdout,'(a,i3,a,1x,f8.3,1x)') 'The eigenvalue number',ii,'is negative and equals to',eigenvalues(ii)
579 0 : write(Invar%stdout,'(a)') 'The Born-Huang stability criterion is not fulfilled'
580 : end if
581 : end do
582 :
583 44 : ABI_FREE(WORK)
584 44 : ABI_FREE(Sij)
585 44 : ABI_FREE(eigenvalues)
586 :
587 : ! For an anisotropic material
588 44 : write(Invar%stdout,'(a)') ' '
589 44 : write(Invar%stdout,'(a)') ' ========== For an Anisotropic Material ======================================='
590 1892 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
591 308 : ABI_MALLOC(IPIV,(6)); IPIV(:)=0
592 308 : ABI_MALLOC(WORK,(6)); WORK(:)=0.d0
593 1892 : Sij(:,:)=Cij(:,:)
594 44 : call DGETRF(6,6,Sij,6,IPIV,INFO)
595 44 : call DGETRI(6,Sij,6,IPIV,WORK,6,INFO)
596 44 : ABI_FREE(IPIV)
597 44 : ABI_FREE(WORK)
598 44 : write(Invar%stdout,'(a)') ' Sijkl [in GPa-1]='
599 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S11 S12 S13 S14 S15 S16 | ',Sij(1,1),Sij(1,2),Sij(1,3),Sij(1,4),Sij(1,5),Sij(1,6)
600 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S21 S22 S23 S24 S25 S26 | ',Sij(2,1),Sij(2,2),Sij(2,3),Sij(2,4),Sij(2,5),Sij(2,6)
601 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S31 S32 S33 S34 S35 S36 | ',Sij(3,1),Sij(3,2),Sij(3,3),Sij(3,4),Sij(3,5),Sij(3,6)
602 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S41 S42 S43 S44 S45 S46 | = ',Sij(4,1),Sij(4,2),Sij(4,3),Sij(4,4),Sij(4,5),Sij(4,6)
603 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S51 S52 S53 S54 S55 S56 | ',Sij(5,1),Sij(5,2),Sij(5,3),Sij(5,4),Sij(5,5),Sij(5,6)
604 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S61 S62 S63 S64 S65 S66 | ',Sij(6,1),Sij(6,2),Sij(6,3),Sij(6,4),Sij(6,5),Sij(6,6)
605 1892 : Lattice%Sij=Sij
606 44 : ABI_FREE(Sij)
607 :
608 : ! For an orthotropic material
609 44 : write(Invar%stdout,'(a)') ' '
610 44 : write(Invar%stdout,'(a)') ' ========== For an Orthotropic Material (see B. M. Lempriere (1968)) =========='
611 :
612 : ! Mean value of the off-diagonal elements
613 44 : Cij(1,2)=(Cij(1,2)+Cij(2,1))/2.d0 ; Cij(2,1)=Cij(1,2)
614 44 : Cij(1,3)=(Cij(1,3)+Cij(3,1))/2.d0 ; Cij(3,1)=Cij(1,3)
615 44 : Cij(2,3)=(Cij(2,3)+Cij(3,2))/2.d0 ; Cij(3,2)=Cij(2,3)
616 :
617 : ! Young's modulus
618 : E1=(Cij(1,1)*Cij(2,2)*Cij(3,3)+2.d0*Cij(2,3)*Cij(1,2)*Cij(1,3)-Cij(1,1)*Cij(2,3)**2-Cij(2,2)*Cij(1,3)**2-Cij(3,3)*Cij(1,2)**2)&
619 44 : & /(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
620 : E2=(Cij(1,1)*Cij(2,2)*Cij(3,3)+2.d0*Cij(2,3)*Cij(1,2)*Cij(1,3)-Cij(1,1)*Cij(2,3)**2-Cij(2,2)*Cij(1,3)**2-Cij(3,3)*Cij(1,2)**2)&
621 44 : & /(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
622 : E3=(Cij(1,1)*Cij(2,2)*Cij(3,3)+2.d0*Cij(2,3)*Cij(1,2)*Cij(1,3)-Cij(1,1)*Cij(2,3)**2-Cij(2,2)*Cij(1,3)**2-Cij(3,3)*Cij(1,2)**2)&
623 44 : & /(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
624 44 : write(Invar%stdout,'(a,3(f8.3,1x))') ' Young modulus E1, E2 and E3 [in GPa]=',E1,E2,E3
625 :
626 : ! Poisson Ratio
627 44 : Nu21=(Cij(1,2)*Cij(3,3)-Cij(1,3)*Cij(2,3))/(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
628 44 : Nu31=(Cij(1,3)*Cij(2,2)-Cij(1,2)*Cij(2,3))/(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
629 44 : Nu23=(Cij(1,1)*Cij(2,3)-Cij(1,2)*Cij(1,3))/(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
630 44 : Nu12=(Cij(1,2)*Cij(3,3)-Cij(1,3)*Cij(2,3))/(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
631 44 : Nu13=(Cij(2,2)*Cij(1,3)-Cij(1,2)*Cij(2,3))/(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
632 44 : Nu32=(Cij(1,1)*Cij(2,3)-Cij(1,2)*Cij(1,3))/(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
633 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' Poisson ratio Nu21, Nu31, Nu23, Nu12, Nu13 and Nu32=',Nu21,Nu31,Nu23,Nu12,Nu13,Nu32
634 :
635 : ! Shear modulus
636 44 : G23=Cij(4,4) ; G13=Cij(5,5) ; G12=Cij(6,6)
637 44 : write(Invar%stdout,'(a,3(f8.3,1x))') ' Shear modulus G23, G13 and G12 [in GPa]=',G23,G13,G12
638 :
639 : ! Compliance matrix
640 1892 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
641 44 : Sij(1,1)= 1.d0/E1 ; Sij(1,2)=-Nu21/E2 ; Sij(1,3)=-Nu31/E3 ; Sij(1,4)=0.d0 ; Sij(1,5)=0.d0 ; Sij(1,6)=0.d0
642 44 : Sij(2,1)=-Nu12/E1 ; Sij(2,2)= 1.d0/E2 ; Sij(2,3)=-Nu32/E3 ; Sij(2,4)=0.d0 ; Sij(2,5)=0.d0 ; Sij(2,6)=0.d0
643 44 : Sij(3,1)=-Nu13/E1 ; Sij(3,2)=-Nu23/E2 ; Sij(3,3)= 1.d0/E3 ; Sij(3,4)=0.d0 ; Sij(3,5)=0.d0 ; Sij(3,6)=0.d0
644 44 : Sij(4,1)= 0.d0 ; Sij(4,2)= 0.d0 ; Sij(4,3)= 0.d0 ; Sij(4,4)=1.d0/G23 ; Sij(4,5)=0.d0 ; Sij(4,6)=0.d0
645 44 : Sij(5,1)= 0.d0 ; Sij(5,2)= 0.d0 ; Sij(5,3)= 0.d0 ; Sij(5,4)=0.d0 ; Sij(5,5)=1.d0/G13 ; Sij(5,6)=0.d0
646 44 : Sij(6,1)= 0.d0 ; Sij(6,2)= 0.d0 ; Sij(6,3)= 0.d0 ; Sij(6,4)=0.d0 ; Sij(6,5)=0.d0 ; Sij(6,6)=1.d0/G12
647 : ! Remove the rounding errors before writing (for non regression testing purposes)
648 308 : do ii=1,6
649 1892 : do jj=1,6
650 1848 : if (abs(Sij(ii,jj)).lt.tol8) Sij(ii,jj)=zero
651 : end do
652 : end do
653 44 : write(Invar%stdout,'(a)') ' Sijkl [in GPa-1]='
654 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S11 S12 S13 S14 S15 S16 | ',Sij(1,1),Sij(1,2),Sij(1,3),Sij(1,4),Sij(1,5),Sij(1,6)
655 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S21 S22 S23 S24 S25 S26 | ',Sij(2,1),Sij(2,2),Sij(2,3),Sij(2,4),Sij(2,5),Sij(2,6)
656 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S31 S32 S33 S34 S35 S36 | ',Sij(3,1),Sij(3,2),Sij(3,3),Sij(3,4),Sij(3,5),Sij(3,6)
657 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S41 S42 S43 S44 S45 S46 | = ',Sij(4,1),Sij(4,2),Sij(4,3),Sij(4,4),Sij(4,5),Sij(4,6)
658 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S51 S52 S53 S54 S55 S56 | ',Sij(5,1),Sij(5,2),Sij(5,3),Sij(5,4),Sij(5,5),Sij(5,6)
659 44 : write(Invar%stdout,'(a,6(f8.3,1x))') ' | S61 S62 S63 S64 S65 S66 | ',Sij(6,1),Sij(6,2),Sij(6,3),Sij(6,4),Sij(6,5),Sij(6,6)
660 :
661 : !==========================================================================================
662 : !===================== Bulk and Shear modulus--Sound velocities ===========================
663 : !==========================================================================================
664 : ! Voigt notation
665 44 : write(Invar%stdout,'(a,f9.3)')' For density rho [in kg.m-3]=',rho
666 44 : write(Invar%stdout,'(a)')' '
667 44 : write(Invar%stdout,'(a)')' ========================= Voigt average (constant strain) ==================='
668 44 : BV=((Cij(1,1)+Cij(2,2)+Cij(3,3))+2.d0*(Cij(1,2)+Cij(1,3)+Cij(2,3)))/9.d0
669 44 : GV=((Cij(1,1)+Cij(2,2)+Cij(3,3))- (Cij(1,2)+Cij(1,3)+Cij(2,3))+3.d0*(Cij(4,4)+Cij(5,5)+Cij(6,6)))/15.d0
670 44 : write(Invar%stdout,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BV,' and Shear G=',GV
671 44 : Eaverage=9.d0*BV*GV/(3*BV+GV)
672 44 : Nuaverage=0.5*(1.d0-(3.d0*GV)/(3.d0*BV+GV) )
673 44 : Laverage=(3.d0*BV-2.d0*GV)/3.d0
674 44 : write(Invar%stdout,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
675 88 : & ' and Poisson ratio Nu=',Nuaverage
676 44 : Vp=dsqrt(1.d9*(BV+4.d0*GV/3.d0)/rho)
677 44 : Vs=dsqrt(1.d9*GV/rho)
678 44 : Vphi=dsqrt(1.d9*BV/rho)
679 44 : write(Invar%stdout,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
680 44 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
681 44 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Lattice%ucvol)**(1./3.)*Ha_K
682 44 : write(Invar%stdout,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
683 :
684 : ! Reuss notation
685 44 : write(Invar%stdout,'(a)')' '
686 44 : write(Invar%stdout,'(a)')' ========================= Reuss average (constant stress) ==================='
687 44 : BR=1.d0/(Sij(1,1)+Sij(2,2)+Sij(3,3)+2.d0*(Sij(1,2)+Sij(1,3)+Sij(2,3)))
688 44 : GR=15.d0/(4.d0*(Sij(1,1)+Sij(2,2)+Sij(3,3))-4.d0*(Sij(1,2)+Sij(1,3)+Sij(2,3))+3.d0*(Sij(4,4)+Sij(5,5)+Sij(6,6)))
689 44 : write(Invar%stdout,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BR,' and Shear G=',GR
690 44 : Eaverage=9.d0*BR*GR/(3*BR+GR)
691 44 : Nuaverage=0.5*(1.d0-(3.d0*GR)/(3.d0*BR+GR) )
692 44 : Laverage=(3.d0*BR-2.d0*GR)/3.d0
693 44 : write(Invar%stdout,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
694 88 : & ' and Poisson ratio Nu=',Nuaverage
695 44 : Vp=dsqrt(1.d9*(BR+4.d0*GR/3.d0)/rho)
696 44 : Vs=dsqrt(1.d9*GR/rho)
697 44 : Vphi=dsqrt(1.d9*BR/rho)
698 44 : write(Invar%stdout,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
699 44 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
700 44 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Lattice%ucvol)**(1./3.)*Ha_K
701 44 : write(Invar%stdout,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
702 :
703 : ! Voigt-Reuss-Hill notation
704 44 : write(Invar%stdout,'(a)')' '
705 44 : write(Invar%stdout,'(a)')' ============================== Hill average ================================='
706 44 : BH=(BR+BV)/2.d0
707 44 : GH=(GR+GV)/2.d0
708 44 : write(Invar%stdout,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BH,' and Shear G=',GH
709 44 : Eaverage=9.d0*BH*GH/(3*BH+GH)
710 44 : Nuaverage=0.5*(1.d0-(3.d0*GH)/(3.d0*BH+GH) )
711 44 : Laverage=(3.d0*BH-2.d0*GH)/3.d0
712 44 : write(Invar%stdout,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
713 88 : & ' and Poisson ratio Nu=',Nuaverage
714 44 : Vp=dsqrt(1.d9*(BH+4.d0*GH/3.d0)/rho)
715 44 : Vs=dsqrt(1.d9*GH/rho)
716 44 : Vphi=dsqrt(1.d9*BH/rho)
717 44 : write(Invar%stdout,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
718 44 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
719 44 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Lattice%ucvol)**(1./3.)*Ha_K
720 44 : write(Invar%stdout,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
721 :
722 : ! Compute the elastic anisotropy
723 44 : A_U=5.d0*GV/GR+BV/BR -6.d0
724 44 : A_B=(BV-BR)/(BV+BR)
725 44 : A_G=(GV-GR)/(GV+GR)
726 44 : write(Invar%stdout,'(a)')' '
727 44 : write(Invar%stdout,'(a)')' ========================= Elastic anisotropy ================================='
728 44 : write(Invar%stdout,'(a,f9.3)')' Elastic anisotropy index : A_U= 5*G_V/G_R + K_V/K_R - 6 =',A_U
729 44 : write(Invar%stdout,'(a,f9.3)')' Bulk anisotropy ratio : A_B= (B_V-B_R)/(B_V+B_R) =',A_B
730 44 : write(Invar%stdout,'(a,f9.3)')' Shear anisotropy ratio : A_G= (G_V-G_R)/(G_V+G_R) =',A_G
731 :
732 :
733 : ! Store the RVH value of the elastic moduli (will be useful for the Gruneisen)
734 44 : Lattice%BulkModulus_T=BH
735 44 : Lattice%Shear=GH
736 44 : Lattice%Density=rho
737 :
738 44 : ABI_FREE(cijkl)
739 44 : ABI_FREE(Cij)
740 44 : ABI_FREE(Sij)
741 44 : end subroutine tdep_calc_elastic
742 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
743 :
744 : end module m_tdep_phdos
|
