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_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
|
