Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_model
2 : : !! NAME
3 : : !! m_tdep_model
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Model data type
7 : : !! which holds the IFC at all orders and derived quantities.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_tdep_model
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : : use m_abicore
28 : : use m_xmpi
29 : : use m_io_tools
30 : : use m_numeric_tools
31 : : use m_matrix, only : mat33det
32 : : use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt
33 : : use m_phonons, only : phdos_t
34 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
35 : : use m_tdep_latt, only : Lattice_type
36 : : use m_tdep_qpt, only : Qpoints_type, Qbz_type
37 : : use m_tdep_sampling, only : tdep_Sampling_type
38 : : use m_tdep_shell, only : Shell_type
39 : : use m_tdep_sym, only : Symmetries_type
40 : : use m_tdep_phi2, only : Phi2_type, Eigen_type, tdep_write_phi1, tdep_write_phi2
41 : : use m_tdep_phi3, only : tdep_write_phi3, tdep_build_phi3_333
42 : : use m_tdep_phi4, only : tdep_build_phi4_3333, tdep_write_phi4
43 : : use m_tdep_supercell, only : tdep_Supercell_type
44 : : use m_tdep_constraints, only : Constraints_type
45 : :
46 : : implicit none
47 : : !!***
48 : :
49 : : !----------------------------------------------------------------------
50 : :
51 : : !!****t* m_tdep_model/tdep_Model_type
52 : : !! NAME
53 : : !! tdep_Model_type
54 : : !!
55 : : !! FUNCTION
56 : : !! Holds the result of the TDEP fit, that is, the IFCs and resulting forces
57 : : !! and energy contributions.
58 : : !!
59 : : !! SOURCE
60 : :
61 : : type,public :: tdep_Model_type
62 : :
63 : : integer :: order
64 : : ! Maximum IFC order (2, 3, or 4)
65 : :
66 : : integer :: natom
67 : : ! Number of atoms
68 : :
69 : : integer :: natom_unitcell
70 : : ! Number of atoms in the unitcell
71 : :
72 : : integer :: my_nstep
73 : : ! Number of MD steps held locally
74 : :
75 : : integer :: nshell3rd
76 : : ! Number of shells for 3rd order IFC
77 : :
78 : : integer :: nshell4th
79 : : ! Number of shells for 4th order IFC
80 : :
81 : : double precision :: U0
82 : : ! Constant term contribution to energy
83 : :
84 : : double precision :: temperature
85 : : ! Temperature of the MD simulation
86 : :
87 : : double precision :: ucvol
88 : : ! Unit cell volume (bohr^3), per unitcell
89 : :
90 : : double precision :: Sij(6,6)
91 : : ! Elastic compliance tensor, in Voigt notation (GPa^-1)
92 : :
93 : : double precision :: BulkModulus_T
94 : : ! Isothermal bulk modulus (GPa)
95 : :
96 : : double precision :: BulkModulus_S
97 : : ! Isentropic (adiabatic) bulk modulus (GPa)
98 : :
99 : : double precision :: HeatCapa_V
100 : : ! Heat capacity at constant volume (k_B/f.u.)
101 : :
102 : : double precision :: HeatCapa_P
103 : : ! Heat capacity at constant pressure (k_B/f.u.)
104 : :
105 : : double precision :: Shear
106 : : ! Hill-averaged shear modulus (GPa)
107 : :
108 : : double precision :: Density
109 : : ! Mass density (kg/m^3)
110 : :
111 : : double precision, allocatable :: Forces(:)
112 : : ! Forces(3*natom*my_nstep)
113 : : ! Forces computed from the fitted IFC and the displacements
114 : :
115 : : double precision, allocatable :: Phi1(:)
116 : : ! Phi1(3*natom)
117 : : ! First-order IFC
118 : :
119 : : type(Phi2_type) :: Phi2
120 : : ! Second-order IFC
121 : :
122 : : double precision, allocatable :: Phi3(:,:,:,:)
123 : : ! Phi3(3,3,3,Shell3at%nshell)
124 : : ! Third-order IFC
125 : :
126 : : double precision, allocatable :: Phi4(:,:,:,:,:)
127 : : ! Phi3(3,3,3,3,Shell4at%nshell)
128 : : ! Fourth-order IFC
129 : :
130 : : double precision, allocatable :: Phi1Ui(:)
131 : : ! Phi1Ui(my_nstep)
132 : : ! First-order IFC contribution to energy
133 : :
134 : : double precision, allocatable :: Phi2UiUj(:)
135 : : ! Phi2UiUj(my_nstep)
136 : : ! Second-order IFC contribution to energy
137 : :
138 : : double precision, allocatable :: Phi3UiUjUk(:)
139 : : ! Phi3UiUjUk(my_nstep)
140 : : ! Third-order IFC contribution to energy
141 : :
142 : : double precision, allocatable :: Phi4UiUjUkUl(:)
143 : : ! Phi4UiUjUkUl(my_nstep)
144 : : ! Fourth-order IFC contribution to energy
145 : :
146 : : contains
147 : :
148 : : procedure :: init => tdep_model_init
149 : : procedure :: free => tdep_model_free
150 : : procedure :: compute_energy => tdep_model_compute_energy
151 : : procedure :: calc_elastic => tdep_model_calc_elastic
152 : : procedure :: calc_thermo => tdep_model_calc_thermo
153 : : procedure :: calc_ftot2 => tdep_calc_ftot2
154 : : procedure :: calc_ftot3 => tdep_calc_ftot3
155 : : procedure :: calc_ftot4 => tdep_calc_ftot4
156 : : procedure :: write_phi1 => tdep_model_write_phi1
157 : : procedure :: write_phi2 => tdep_model_write_phi2
158 : : procedure :: write_phi3 => tdep_model_write_phi3
159 : : procedure :: write_phi4 => tdep_model_write_phi4
160 : : procedure :: check_constraints => tdep_check_constraints
161 : : procedure :: calc_alpha_gamma => tdep_calc_alpha_gamma
162 : : procedure :: calc_gruneisen => tdep_calc_write_gruneisen
163 : :
164 : : end type tdep_Model_type
165 : : !!***
166 : :
167 : : public :: tdep_calc_ftot2
168 : : public :: tdep_calc_ftot3
169 : : public :: tdep_calc_ftot4
170 : : public :: tdep_model_write_phi1
171 : : public :: tdep_model_write_phi2
172 : : public :: tdep_model_write_phi3
173 : : public :: tdep_model_write_phi4
174 : : public :: tdep_check_constraints
175 : : public :: tdep_calc_gruneisen
176 : : public :: tdep_calc_write_gruneisen
177 : : public :: tdep_calc_alpha_gamma
178 : :
179 : : contains
180 : :
181 : : !=====================================================================================================
182 : :
183 : : !!****f* m_tdep_model/tdep_model_init
184 : : !! NAME
185 : : !! tdep_model_init
186 : : !!
187 : : !! FUNCTION
188 : : !! Initialize the TDEP Model object: set the dimensions, compute the unit
189 : : !! cell volume, and allocate (and set to zero) the IFC at all orders
190 : : !! as well as the arrays holding the forces and energy contributions.
191 : : !!
192 : : !! INPUTS
193 : : !! Invar = Input object containing the input variables.
194 : : !! Lattice = Lattice object describing the ideal structure.
195 : : !! Shell3at = Shell object for third-order interactions
196 : : !! (only used if Invar%order >= 3).
197 : : !! Shell4at = Shell object for fourth-order interactions
198 : : !! (only used if Invar%order == 4).
199 : : !!
200 : : !! OUTPUT
201 : : !! Model = TDEP Model object, initialized.
202 : : !!
203 : : !! SOURCE
204 : :
205 [ + - ]: 46 : subroutine tdep_model_init(Model, Invar, Lattice, Shell3at, Shell4at)
206 : :
207 : : class(tdep_Model_type),intent(out) :: Model
208 : : type(atdep_dataset_type),intent(in) :: Invar
209 : : type(Lattice_type),intent(in) :: Lattice
210 : : type(Shell_type),intent(in) :: Shell3at, Shell4at
211 : :
212 : 46 : Model%order = Invar%order
213 : 46 : Model%natom = Invar%natom
214 : 46 : Model%natom_unitcell = Invar%natom_unitcell
215 : 46 : Model%my_nstep = Invar%my_nstep
216 : 46 : Model%temperature = Invar%temperature
217 : 46 : Model%nshell3rd = 1
218 : 46 : Model%nshell4th = 1
219 : :
220 : 46 : Model%U0 = zero
221 : :
222 : : ! Unit cell volume: computed here since Lattice is available; the other
223 : : ! quantities derived from the model (elastic constants, thermodynamics)
224 : : ! are computed and stored later on, once they become available.
225 : 46 : Model%ucvol = abs(mat33det(Lattice%rprimd))
226 [ + + + + ]: 1978 : Model%Sij(:,:) = zero
227 : 46 : Model%BulkModulus_T = zero
228 : 46 : Model%BulkModulus_S = zero
229 : 46 : Model%HeatCapa_V = zero
230 : 46 : Model%HeatCapa_P = zero
231 : 46 : Model%Shear = zero
232 : 46 : Model%Density = zero
233 : :
234 [ + - - + : 441354 : ABI_CALLOC(Model%Forces, (3*Model%natom*Model%my_nstep))
+ - + - +
- - + - +
- + + + ]
235 [ - + + - : 1130 : ABI_CALLOC(Model%Phi1Ui ,(Model%my_nstep))
+ - + - -
+ - + - +
+ + ]
236 [ - + + - : 1084 : ABI_CALLOC(Model%Phi2UiUj ,(Model%my_nstep))
- + - + -
+ + + ]
237 [ - + + - : 1084 : ABI_CALLOC(Model%Phi3UiUjUk ,(Model%my_nstep))
- + - + -
+ + + ]
238 [ - + + - : 1084 : ABI_CALLOC(Model%Phi4UiUjUkUl,(Model%my_nstep))
- + - + -
+ + + ]
239 : :
240 [ + - - + : 21972 : ABI_CALLOC(Model%Phi1,(3*Model%natom))
+ - + - +
- - + - +
- + + + ]
241 : :
242 : 46 : call Model%Phi2%init(Invar%dipdip,Invar%natom)
243 : :
244 [ + + ]: 46 : if (Model%order.ge.3) then
245 : 16 : Model%nshell3rd = Shell3at%nshell
246 [ - + + - : 2328 : ABI_CALLOC(Model%Phi3,(3,3,3,Model%nshell3rd))
+ - + - -
+ - + - +
+ + + + +
+ + + ]
247 : : end if
248 : :
249 [ + + ]: 46 : if (Model%order.eq.4) then
250 : 9 : Model%nshell4th = Shell4at%nshell
251 [ - + + - : 5956 : ABI_CALLOC(Model%Phi4,(3,3,3,3,Model%nshell4th))
+ - + - -
+ - + - +
+ + + + +
+ + + +
+ ]
252 : : end if
253 : :
254 : 46 : end subroutine tdep_model_init
255 : : !!***
256 : :
257 : : !=====================================================================================================
258 : :
259 : : !!****f* m_tdep_model/tdep_model_free
260 : : !! NAME
261 : : !! tdep_model_free
262 : : !!
263 : : !! FUNCTION
264 : : !! Free the memory allocated in the TDEP Model object.
265 : : !!
266 : : !! SIDE EFFECTS
267 : : !! Model = TDEP Model object. All its allocated arrays are deallocated.
268 : : !!
269 : : !! SOURCE
270 : :
271 : 46 : subroutine tdep_model_free(Model)
272 : :
273 : : class(tdep_Model_type),intent(inout) :: Model
274 : :
275 [ - + ]: 46 : ABI_FREE(Model%Forces)
276 [ - + ]: 46 : ABI_FREE(Model%Phi1Ui)
277 [ - + ]: 46 : ABI_FREE(Model%Phi2UiUj)
278 [ - + ]: 46 : ABI_FREE(Model%Phi3UiUjUk)
279 [ - + ]: 46 : ABI_FREE(Model%Phi4UiUjUkUl)
280 : :
281 [ - + ]: 46 : ABI_FREE(Model%Phi1)
282 : :
283 : 46 : call Model%Phi2%free()
284 : :
285 [ + + ]: 46 : if (Model%order.ge.3) then
286 [ - + ]: 16 : ABI_FREE(Model%Phi3)
287 : : end if
288 : :
289 [ + + ]: 46 : if (Model%order.eq.4) then
290 [ - + ]: 9 : ABI_FREE(Model%Phi4)
291 : : end if
292 : :
293 : 46 : end subroutine tdep_model_free
294 : : !!***
295 : :
296 : : !====================================================================================================
297 : :
298 : : !!****f* m_tdep_model/tdep_model_compute_energy
299 : : !! NAME
300 : : !! tdep_model_compute_energy
301 : : !!
302 : : !! FUNCTION
303 : : !! Compute the different contributions to the total energy of the model
304 : : !! (U_0, U_1, U_2, U_3, U_4), averaged over the MD trajectory, together with
305 : : !! the errors of the model with respect to the MD energies and forces
306 : : !! (Delta_U, Delta_U2, Delta_F2, Sigma). Write them in the output file.
307 : : !! Also write the energies (U_MD vs U_TDEP) in the _etotMDvsTDEP.dat file
308 : : !! and the forces (F_MD vs F_TDEP) in the _fcartMDvsTDEP.dat file.
309 : : !!
310 : : !! INPUTS
311 : : !! Invar = Input object containing the input variables.
312 : : !! MPIdata = Info on MPI parallelism.
313 : : !! MD = TDEP Sampling object containing the MD energies, forces and weights.
314 : : !!
315 : : !! OUTPUT
316 : : !!
317 : : !! SIDE EFFECTS
318 : : !! Model = TDEP Model object. On entry, Model%Forces and the energy
319 : : !! contributions (Model%Phi1Ui, Model%Phi2UiUj, Model%Phi3UiUjUk and
320 : : !! Model%Phi4UiUjUkUl) must already be computed (see tdep_calc_ftot2,
321 : : !! tdep_calc_ftot3 and tdep_calc_ftot4). On exit, Model%U0 is
322 : : !! computed and stored.
323 : : !!
324 : : !! SOURCE
325 : :
326 : 46 : subroutine tdep_model_compute_energy(Model,Invar,MPIdata,MD)
327 : :
328 : : class(tdep_Model_type),intent(inout) :: Model
329 : : type(atdep_dataset_type),intent(in) :: Invar
330 : : type(MPI_enreg_type), intent(in) :: MPIdata
331 : : type(tdep_Sampling_type),intent(in) :: MD
332 : :
333 : : integer :: ii,jj,istep,iatom
334 : : double precision :: Delta_F2,Delta_U,Delta_U2
335 : : double precision :: sigma,U_1,U_2,U_3,U_4,UMD
336 : 46 : double precision, allocatable :: tmp(:),Phi_tot(:)
337 : 46 : double precision, allocatable :: U_MD(:),U_TDEP(:),weights_tot(:)
338 : : integer :: ierr
339 : :
340 : : ! Compute the different contributions to total energy from the model
341 : :
342 : 46 : write(ab_out,*)' '
343 : 46 : write(ab_out,*) '#############################################################################'
344 : 46 : write(ab_out,*) '######################### Energies, errors,... #############################'
345 : 46 : write(ab_out,*) '#############################################################################'
346 : :
347 : : ! Compute U0, U_TDEP, Delta_U and write them in the output file
348 : 46 : write(ab_out,'(a)') ' Thermodynamic quantities and convergence parameters of THE MODEL,'
349 : 46 : write(ab_out,'(a)') ' as a function of the step number (energies in eV/atom and forces in Ha/bohr) :'
350 [ + + ]: 46 : if (Model%order.eq.4) then
351 : 9 : write(ab_out,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2 + U_3 + U_4'
352 : 9 : write(ab_out,'(2a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj ',&
353 : 18 : & '- 1/6 sum_ijk Phi3 ui uj uk - 1/24 sum_ijkl Phi4 ui uj uk ul >'
354 : 9 : write(ab_out,'(a)') ' and U_1 = < sum_i Phi1 ui >'
355 : 9 : write(ab_out,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
356 : 9 : write(ab_out,'(a)') ' and U_3 = < 1/6 sum_ijk Phi3 ui uj uk >'
357 : 9 : write(ab_out,'(a)') ' and U_4 = < 1/24 sum_ijkl Phi4 ui uj uk ul >'
358 [ + + ]: 37 : else if (Model%order.eq.3) then
359 : 7 : write(ab_out,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2 + U_3'
360 : 7 : write(ab_out,'(a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj - 1/6 sum_ijk Phi3 ui uj uk >'
361 : 7 : write(ab_out,'(a)') ' and U_1 = < sum_i Phi1 ui >'
362 : 7 : write(ab_out,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
363 : 7 : write(ab_out,'(a)') ' and U_3 = < 1/6 sum_ijk Phi3 ui uj uk >'
364 : : else
365 : 30 : write(ab_out,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2'
366 : 30 : write(ab_out,'(a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj >'
367 : 30 : write(ab_out,'(a)') ' and U_1 = < sum_i Phi1 ui >'
368 : 30 : write(ab_out,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
369 : : end if
370 : 46 : write(ab_out,'(a)') ' Delta_U = < U_MD - U_TDEP > '
371 : 46 : write(ab_out,'(a)') ' Delta_U2= (< (U_MD - U_TDEP)^2 >)**0.5 '
372 : 46 : write(ab_out,'(a)') ' Delta_F2= (< (F_MD - F_TDEP)^2 >)**0.5 '
373 : 46 : write(ab_out,'(a)') ' Sigma = (< (F_MD - F_TDEP)^2 >/<F_MD**2>)**0.5 '
374 [ + + ]: 46 : if (Model%order.eq.4) then
375 : 9 : write(ab_out,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
376 : 18 : & ' U_3 U_4 Delta_U Delta_U2 Delta_F2 Sigma'
377 [ + + ]: 37 : else if (Model%order.eq.3) then
378 : 7 : write(ab_out,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
379 : 14 : & ' U_3 Delta_U Delta_U2 Delta_F2 Sigma'
380 : : else
381 : 30 : write(ab_out,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
382 : 60 : & ' Delta_U Delta_U2 Delta_F2 Sigma'
383 : : end if
384 : :
385 : : ! Compute eucledian distance for forces
386 [ - + + + ]: 552 : ABI_MALLOC(tmp,(11)) ; tmp(:) =0.d0
387 [ + + ]: 1038 : do istep=1,MD%my_nstep
388 [ + + ]: 148110 : do iatom=1,MD%natom
389 [ + + ]: 589280 : do ii=1,3
390 : 441216 : jj = ii + 3*(iatom-1) + 3*MD%natom*(istep-1)
391 : 441216 : tmp(4)=tmp(4)+(MD%Forces(jj)-Model%Forces(jj))**2*MD%weights(istep)
392 : 588288 : tmp(5)=tmp(5)+MD%Forces(jj)**2*MD%weights(istep)
393 : : end do
394 : : end do
395 : : end do
396 : : ! Compute energies
397 [ - + + - : 1130 : ABI_MALLOC(U_TDEP, (MD%nstep_tot)) ; U_TDEP(:)=0.d0
+ - + - -
+ - + +
+ ]
398 [ - + + - : 1084 : ABI_MALLOC(U_MD, (MD%nstep_tot)) ; U_MD(:) =0.d0
- + - + +
+ ]
399 [ - + + - : 1084 : ABI_MALLOC(weights_tot,(MD%nstep_tot)) ; weights_tot(:)=0.d0
- + - + +
+ ]
400 [ - + + - : 1130 : ABI_MALLOC(Phi_tot, (MPIdata%my_nstep)); Phi_tot(:)=0.d0
+ - + - -
+ - + +
+ ]
401 [ + + ]: 1038 : do istep=1,MD%my_nstep
402 : 992 : tmp(7) =tmp(7) +MD%etot(istep)*MD%weights(istep)
403 : 992 : tmp(10)=tmp(10)+Model%Phi1Ui(istep)*MD%weights(istep)
404 : 992 : tmp(6) =tmp(6) +Model%Phi2UiUj(istep)*MD%weights(istep)
405 : 992 : tmp(8) =tmp(8) +Model%Phi3UiUjUk(istep)*MD%weights(istep)
406 : 1038 : tmp(11)=tmp(11)+Model%Phi4UiUjUkUl(istep)*MD%weights(istep)
407 : : end do
408 : 46 : call xmpi_sum(tmp,MPIdata%comm_step,ierr)
409 : 46 : tmp(1) = tmp(7)-tmp(10)-tmp(6)-tmp(8)-tmp(11)
410 [ + + ]: 1038 : Phi_tot(:)=tmp(1)+Model%Phi1Ui(:)+Model%Phi2UiUj(:)+Model%Phi3UiUjUk(:)+Model%Phi4UiUjUkUl(:)
411 : : call xmpi_gatherv(Phi_tot,MD%my_nstep,U_TDEP,MPIdata%nstep_all,MPIdata%shft_step,&
412 : 46 : & MPIdata%master,MPIdata%comm_step,ierr)
413 : : call xmpi_gatherv(MD%etot,MD%my_nstep,U_MD,MPIdata%nstep_all,MPIdata%shft_step,&
414 : 46 : & MPIdata%master,MPIdata%comm_step,ierr)
415 : : call xmpi_gatherv(MD%weights,MD%my_nstep,weights_tot,MPIdata%nstep_all,MPIdata%shft_step,&
416 : 46 : & MPIdata%master,MPIdata%comm_step,ierr)
417 [ + + ]: 1038 : do istep=1,MD%nstep_tot
418 : 992 : tmp(2) =tmp(2) + (U_MD(istep)-U_TDEP(istep)) * weights_tot(istep)
419 : 1038 : tmp(9) =tmp(9) + (U_MD(istep)-U_TDEP(istep))**2 * weights_tot(istep)
420 : : end do
421 : 46 : Model%U0 =tmp(1) /real(MD%natom)
422 : 46 : UMD =tmp(7) /real(MD%natom)
423 : 46 : U_1 =tmp(10)/real(MD%natom)
424 : 46 : U_2 =tmp(6) /real(MD%natom)
425 : 46 : U_3 =tmp(8) /real(MD%natom)
426 : 46 : U_4 =tmp(11)/real(MD%natom)
427 : 46 : Delta_U =tmp(2) /real(MD%natom)
428 : 46 : Delta_U2 =tmp(9) /real(MD%natom)
429 : 46 : Delta_F2 =tmp(4) /real(MD%natom*3)
430 [ + + ]: 46 : if (tmp(5).eq.0.d0) then
431 : 1 : sigma = 0.d0
432 : : else
433 : 45 : sigma = dsqrt(tmp(4)/tmp(5))
434 : : end if
435 [ + + ]: 46 : if (Model%order.eq.4) then
436 : 9 : write(ab_out,'(10(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,U_3*Ha_eV,U_4*Ha_eV,&
437 : 18 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
438 [ + + ]: 37 : else if (Model%order.eq.3) then
439 : 7 : write(ab_out,'(9(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,U_3*Ha_eV,&
440 : 14 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
441 : : else
442 : 30 : write(ab_out,'(8(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,&
443 : 60 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
444 : : endif
445 : 46 : ABI_FREE(tmp)
446 : 46 : write(ab_out,'(a,1x,f12.5)') ' NOTE : in the harmonic and classical limit (T>>T_Debye), U_2=3/2*kB*T=',&
447 : 92 : & 3.d0/2.d0*kb_HaK*Ha_eV*Model%temperature
448 : :
449 : : ! Write : i) (U_TDEP vs U_MD) in etotMDvsTDEP.dat
450 : : ! ii) (Model%Forces vs MD%Forces) in fcartMDvsTDEP.dat
451 : 46 : write(ab_out,'(a)') ' '
452 : 46 : write(ab_out,'(a)') ' Writing model total energy in file etotMDvsTDEP.dat'
453 : 46 : write(ab_out,'(a)') ' Writing model forces in file fcartMDvsTDEP.dat'
454 [ + - ]: 46 : if (MPIdata%iam_master) then
455 : 46 : open(unit=32,file=trim(Invar%output_prefix)//'_etotMDvsTDEP.dat')
456 : 46 : open(unit=33,file=trim(Invar%output_prefix)//'_fcartMDvsTDEP.dat')
457 : 46 : write(32,'(a)') '# Istep U_MD(Ha) U_TDEP(Ha)'
458 : 46 : write(33,'(a)') '# Forces_MD(Ha/bohr) Forces_TDEP(Ha/bohr)'
459 [ + + ]: 1038 : do istep=1,MD%nstep_tot
460 : 1038 : write(32,'(i6,1x,2(f17.6,1x))') istep,U_MD(istep),U_TDEP(istep)
461 : : end do
462 [ + + ]: 1038 : do istep=1,MD%my_nstep
463 [ + + ]: 148110 : do iatom=1,MD%natom
464 [ + + ]: 589280 : do ii=1,3
465 : 441216 : write(33,'(2(f17.10,1x))') MD%Forces (ii+3*(iatom-1)+3*MD%natom*(istep-1)),&
466 : 1029504 : & Model%Forces(ii+3*(iatom-1)+3*MD%natom*(istep-1))
467 : : end do
468 : : end do
469 : : end do
470 : 46 : close(32)
471 : 46 : close(33)
472 : : end if
473 [ - + ]: 46 : ABI_FREE(U_MD)
474 [ - + ]: 46 : ABI_FREE(U_TDEP)
475 [ - + ]: 46 : ABI_FREE(Phi_tot)
476 [ - + ]: 46 : ABI_FREE(weights_tot)
477 : :
478 : 92 : end subroutine tdep_model_compute_energy
479 : : !!***
480 : :
481 : : !====================================================================================================
482 : :
483 : : !!****f* m_tdep_model/tdep_model_calc_elastic
484 : : !! NAME
485 : : !! tdep_model_calc_elastic
486 : : !!
487 : : !! FUNCTION
488 : : !! Compute the elastic constants (Cijkl and their inverse, the compliance
489 : : !! tensor Sij), the density, and the Voigt/Reuss/Hill-averaged bulk and
490 : : !! shear moduli and sound velocities, from the harmonic (2nd order) IFC.
491 : : !! See Wallace, "Statistical physics of crystals and liquids, World
492 : : !! Scientific" (eqs. 12.28 and 12.29) and B. M. Lempriere (1968) for the
493 : : !! orthotropic-material formulas.
494 : : !!
495 : : !! INPUTS
496 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
497 : : !! Invar = Input object containing the input variables.
498 : : !!
499 : : !! OUTPUT
500 : : !!
501 : : !! SIDE EFFECTS
502 : : !! Model = TDEP Model object. On entry, Model%Phi2%SR and Model%ucvol must
503 : : !! already be set. On exit, Model%Sij, Model%BulkModulus_T,
504 : : !! Model%Shear and Model%Density are computed and stored.
505 : : !!
506 : : !! SOURCE
507 : :
508 : 46 : subroutine tdep_model_calc_elastic(Model,SC,Invar)
509 : :
510 : : integer :: iatom,ii,jj,kk,ll,iatcell,itypat
511 : : integer :: INFO,LWORK
512 : : double precision :: BH,BR,BV,GR,GV,GH,Eaverage,Nuaverage,Laverage,Vp,Vs,Vphi
513 : : double precision :: rho,E1,E2,E3,Nu12,Nu13,Nu23,Nu21,Nu31,Nu32,G23,G13,G12
514 : : double precision :: mass_amu,bohr,V_Debye,T_Debye,A_U,A_B,A_G
515 : 46 : integer, allocatable :: IPIV(:)
516 : 46 : double precision, allocatable :: eigenvalues(:)
517 : 46 : double precision, allocatable :: WORK(:)
518 : 46 : double precision, allocatable :: Sij(:,:),Cij(:,:),aijkl(:,:,:,:),cijkl(:,:,:,:)
519 : : class(tdep_Model_type), intent(inout) :: Model
520 : : type(atdep_dataset_type), intent(in) :: Invar
521 : : type(tdep_Supercell_type), intent(in) :: SC
522 : :
523 : 46 : write(ab_out,'(a)')' '
524 : 46 : write(ab_out,'(a)') ' #############################################################################'
525 : 46 : write(ab_out,'(a)') ' ######################### Elastic constants #################################'
526 : 46 : write(ab_out,'(a)') ' ################ Bulk and Shear modulus--Sound velocities ###################'
527 : 46 : write(ab_out,'(a)') ' #############################################################################'
528 : :
529 : 46 : bohr=Bohr_Ang*1.0e-10
530 : : ! Define atomic mass average
531 : 46 : mass_amu=zero
532 [ + + ]: 156 : do iatom=1,Invar%natom_unitcell
533 : 110 : itypat=Invar%typat_unitcell(iatom)
534 : 156 : mass_amu=mass_amu+Invar%amu(itypat)
535 : : end do
536 : 46 : mass_amu=mass_amu/real(Invar%natom_unitcell)
537 : :
538 : 46 : rho=(mass_amu/1e3)*Invar%natom_unitcell/Model%ucvol/bohr**3/Avogadro
539 : :
540 : : !==========================================================================================
541 : : !===================== Elastic constants ==================================================
542 : : !==========================================================================================
543 : : ! New calculation of elastic constants using the formula (12.28 and 12.29 of
544 : : ! Wallace, Statistical physics of crystals and liquids, Worl Scientific)
545 [ - + + + : 5566 : ABI_MALLOC(aijkl,(3,3,3,3)); aijkl(:,:,:,:)=0.d0
+ + + + +
+ ]
546 [ - + + + : 5566 : ABI_MALLOC(cijkl,(3,3,3,3)); cijkl(:,:,:,:)=0.d0
+ + + + +
+ ]
547 [ + + ]: 184 : do ii=1,3
548 [ + + ]: 598 : do jj=1,3
549 [ + + ]: 1794 : do kk=1,3
550 [ + + ]: 5382 : do ll=1,3
551 [ + + ]: 594486 : do iatom=1,Invar%natom
552 [ + + ]: 1863162 : do iatcell=1,Invar%natom_unitcell
553 : : aijkl(ii,jj,kk,ll)=aijkl(ii,jj,kk,ll)-Model%Phi2%SR(ii+(iatcell-1)*3,3*(iatom-1)+jj)&
554 : 1859436 : & *SC%dij(iatcell,iatom,kk+1)*SC%dij(iatcell,iatom,ll+1)/2.d0/Model%ucvol
555 : : end do
556 : : end do
557 : : enddo
558 : : enddo
559 : : enddo
560 : : enddo
561 : :
562 [ + + ]: 184 : do ii=1,3
563 [ + + ]: 598 : do jj=1,3
564 [ + + ]: 1794 : do kk=1,3
565 [ + + ]: 5382 : do ll=1,3
566 : 4968 : cijkl(ii,jj,kk,ll)=aijkl(ii,kk,jj,ll)+aijkl(jj,kk,ii,ll)-aijkl(ii,jj,kk,ll)
567 : : !FB cijkl(ii,jj,kk,ll)=aijkl(ii,kk,jj,ll)+aijkl(jj,kk,ii,ll)-aijkl(jj,ii,kk,ll)
568 : : enddo
569 : : enddo
570 : : enddo
571 : : enddo
572 [ - + ]: 46 : ABI_FREE(aijkl)
573 : :
574 [ + + + + : 5566 : cijkl(:,:,:,:)=cijkl(:,:,:,:)*29421.033d0
+ + + + ]
575 : :
576 [ - + + + : 1978 : ABI_MALLOC(Cij,(6,6)) ; Cij(:,:)=0.d0
+ + ]
577 : 46 : Cij(1,1)=cijkl(1,1,1,1) ; Cij(1,2)=cijkl(1,1,2,2) ; Cij(1,3)=cijkl(1,1,3,3)
578 : 46 : Cij(1,4)=cijkl(1,1,2,3) ; Cij(1,5)=cijkl(1,1,1,3) ; Cij(1,6)=cijkl(1,1,1,2)
579 : 46 : Cij(2,1)=cijkl(2,2,1,1) ; Cij(2,2)=cijkl(2,2,2,2) ; Cij(2,3)=cijkl(2,2,3,3)
580 : 46 : Cij(2,4)=cijkl(2,2,2,3) ; Cij(2,5)=cijkl(2,2,1,3) ; Cij(2,6)=cijkl(2,2,1,2)
581 : 46 : Cij(3,1)=cijkl(3,3,1,1) ; Cij(3,2)=cijkl(3,3,2,2) ; Cij(3,3)=cijkl(3,3,3,3)
582 : 46 : Cij(3,4)=cijkl(3,3,2,3) ; Cij(3,5)=cijkl(3,3,1,3) ; Cij(3,6)=cijkl(3,3,1,2)
583 : 46 : Cij(4,1)=cijkl(2,3,1,1) ; Cij(4,2)=cijkl(2,3,2,2) ; Cij(4,3)=cijkl(2,3,3,3)
584 : 46 : Cij(4,4)=cijkl(2,3,2,3) ; Cij(4,5)=cijkl(2,3,1,3) ; Cij(4,6)=cijkl(2,3,1,2)
585 : 46 : Cij(5,1)=cijkl(1,3,1,1) ; Cij(5,2)=cijkl(1,3,2,2) ; Cij(5,3)=cijkl(1,3,3,3)
586 : 46 : Cij(5,4)=cijkl(1,3,2,3) ; Cij(5,5)=cijkl(1,3,1,3) ; Cij(5,6)=cijkl(1,3,1,2)
587 : 46 : Cij(6,1)=cijkl(1,2,1,1) ; Cij(6,2)=cijkl(1,2,2,2) ; Cij(6,3)=cijkl(1,2,3,3)
588 : 46 : Cij(6,4)=cijkl(1,2,2,3) ; Cij(6,5)=cijkl(1,2,1,3) ; Cij(6,6)=cijkl(1,2,1,2)
589 : : ! Remove the rounding errors before writing (for non regression testing purposes)
590 [ + + ]: 322 : do ii=1,6
591 [ + + ]: 1978 : do jj=1,6
592 [ + + ]: 1932 : if (abs(Cij(ii,jj)).lt.tol8) Cij(ii,jj)=zero
593 : : end do
594 : : end do
595 : 46 : write(ab_out,'(a)') ' '
596 : 46 : write(ab_out,'(a)') ' ========== Using the formulation proposed by Wallace (using the IFC) ========='
597 : 46 : write(ab_out,'(a)') ' Cijkl [in GPa]='
598 : 46 : write(ab_out,'(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)
599 : 46 : write(ab_out,'(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)
600 : 46 : write(ab_out,'(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)
601 : 46 : write(ab_out,'(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)
602 : 46 : write(ab_out,'(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)
603 : 46 : write(ab_out,'(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)
604 : :
605 : : ! Compute the eigenvalues of the Cij matrix in order to find the Born-Huang
606 : : ! stability criterion (see Wallace, Thermodynamics of crystal, p39)
607 [ - + ]: 46 : ABI_MALLOC(WORK,(1))
608 [ - + + + : 1978 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
+ + ]
609 [ - + + + ]: 322 : ABI_MALLOC(eigenvalues,(6)) ; eigenvalues(:)=0.d0
610 [ + + + + ]: 1978 : Sij(:,:)=Cij(:,:)
611 : 46 : LWORK=-1
612 : 46 : call DSYEV('N','U',6,Sij,6,eigenvalues,WORK,LWORK,INFO)
613 : 46 : LWORK=WORK(1)
614 : 46 : ABI_FREE(WORK)
615 [ - + + - : 138 : ABI_MALLOC(WORK,(LWORK))
+ - + - -
+ - + ]
616 : 46 : call DSYEV('N','U',6,Sij,6,eigenvalues,WORK,LWORK,INFO)
617 [ + + ]: 322 : do ii=1,6
618 [ - + ]: 322 : if (eigenvalues(ii).lt.0.d0) then
619 : 0 : write(ab_out,'(a)') ' WARNING :'
620 : 0 : write(ab_out,'(a,i3,a,1x,f8.3,1x)') 'The eigenvalue number',ii,'is negative and equals to',eigenvalues(ii)
621 : 0 : write(ab_out,'(a)') 'The Born-Huang stability criterion is not fulfilled'
622 : : end if
623 : : end do
624 : :
625 [ - + ]: 46 : ABI_FREE(WORK)
626 [ - + ]: 46 : ABI_FREE(Sij)
627 [ - + ]: 46 : ABI_FREE(eigenvalues)
628 : :
629 : : ! For an anisotropic material
630 : 46 : write(ab_out,'(a)') ' '
631 : 46 : write(ab_out,'(a)') ' ========== For an Anisotropic Material ======================================='
632 [ - + + + : 1978 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
+ + ]
633 [ - + + + ]: 322 : ABI_MALLOC(IPIV,(6)); IPIV(:)=0
634 [ - + + + ]: 322 : ABI_MALLOC(WORK,(6)); WORK(:)=0.d0
635 [ + + + + ]: 1978 : Sij(:,:)=Cij(:,:)
636 : 46 : call DGETRF(6,6,Sij,6,IPIV,INFO)
637 : 46 : call DGETRI(6,Sij,6,IPIV,WORK,6,INFO)
638 [ - + ]: 46 : ABI_FREE(IPIV)
639 [ - + ]: 46 : ABI_FREE(WORK)
640 : 46 : write(ab_out,'(a)') ' Sijkl [in GPa-1]='
641 : 46 : write(ab_out,'(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)
642 : 46 : write(ab_out,'(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)
643 : 46 : write(ab_out,'(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)
644 : 46 : write(ab_out,'(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)
645 : 46 : write(ab_out,'(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)
646 : 46 : write(ab_out,'(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)
647 [ + + + + ]: 1978 : Model%Sij=Sij
648 [ - + ]: 46 : ABI_FREE(Sij)
649 : :
650 : : ! For an orthotropic material
651 : 46 : write(ab_out,'(a)') ' '
652 : 46 : write(ab_out,'(a)') ' ========== For an Orthotropic Material (see B. M. Lempriere (1968)) =========='
653 : :
654 : : ! Mean value of the off-diagonal elements
655 : 46 : Cij(1,2)=(Cij(1,2)+Cij(2,1))/2.d0 ; Cij(2,1)=Cij(1,2)
656 : 46 : Cij(1,3)=(Cij(1,3)+Cij(3,1))/2.d0 ; Cij(3,1)=Cij(1,3)
657 : 46 : Cij(2,3)=(Cij(2,3)+Cij(3,2))/2.d0 ; Cij(3,2)=Cij(2,3)
658 : :
659 : : ! Young's modulus
660 : : 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)&
661 : 46 : & /(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
662 : : 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)&
663 : 46 : & /(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
664 : : 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)&
665 : 46 : & /(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
666 : 46 : write(ab_out,'(a,3(f8.3,1x))') ' Young modulus E1, E2 and E3 [in GPa]=',E1,E2,E3
667 : :
668 : : ! Poisson Ratio
669 : 46 : Nu21=(Cij(1,2)*Cij(3,3)-Cij(1,3)*Cij(2,3))/(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
670 : 46 : Nu31=(Cij(1,3)*Cij(2,2)-Cij(1,2)*Cij(2,3))/(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
671 : 46 : Nu23=(Cij(1,1)*Cij(2,3)-Cij(1,2)*Cij(1,3))/(Cij(1,1)*Cij(3,3)-Cij(1,3)**2)
672 : 46 : Nu12=(Cij(1,2)*Cij(3,3)-Cij(1,3)*Cij(2,3))/(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
673 : 46 : Nu13=(Cij(2,2)*Cij(1,3)-Cij(1,2)*Cij(2,3))/(Cij(2,2)*Cij(3,3)-Cij(2,3)**2)
674 : 46 : Nu32=(Cij(1,1)*Cij(2,3)-Cij(1,2)*Cij(1,3))/(Cij(1,1)*Cij(2,2)-Cij(1,2)**2)
675 : 46 : write(ab_out,'(a,6(f8.3,1x))') ' Poisson ratio Nu21, Nu31, Nu23, Nu12, Nu13 and Nu32=',Nu21,Nu31,Nu23,Nu12,Nu13,Nu32
676 : :
677 : : ! Shear modulus
678 : 46 : G23=Cij(4,4) ; G13=Cij(5,5) ; G12=Cij(6,6)
679 : 46 : write(ab_out,'(a,3(f8.3,1x))') ' Shear modulus G23, G13 and G12 [in GPa]=',G23,G13,G12
680 : :
681 : : ! Compliance matrix
682 [ - + + + : 1978 : ABI_MALLOC(Sij,(6,6)) ; Sij(:,:)=0.d0
+ + ]
683 : 46 : 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
684 : 46 : 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
685 : 46 : 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
686 : 46 : 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
687 : 46 : 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
688 : 46 : 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
689 : : ! Remove the rounding errors before writing (for non regression testing purposes)
690 [ + + ]: 322 : do ii=1,6
691 [ + + ]: 1978 : do jj=1,6
692 [ + + ]: 1932 : if (abs(Sij(ii,jj)).lt.tol8) Sij(ii,jj)=zero
693 : : end do
694 : : end do
695 : 46 : write(ab_out,'(a)') ' Sijkl [in GPa-1]='
696 : 46 : write(ab_out,'(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)
697 : 46 : write(ab_out,'(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)
698 : 46 : write(ab_out,'(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)
699 : 46 : write(ab_out,'(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)
700 : 46 : write(ab_out,'(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)
701 : 46 : write(ab_out,'(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)
702 : :
703 : : !==========================================================================================
704 : : !===================== Bulk and Shear modulus--Sound velocities ===========================
705 : : !==========================================================================================
706 : : ! Voigt notation
707 : 46 : write(ab_out,'(a,f9.3)')' For density rho [in kg.m-3]=',rho
708 : 46 : write(ab_out,'(a)')' '
709 : 46 : write(ab_out,'(a)')' ========================= Voigt average (constant strain) ==================='
710 : 46 : BV=((Cij(1,1)+Cij(2,2)+Cij(3,3))+2.d0*(Cij(1,2)+Cij(1,3)+Cij(2,3)))/9.d0
711 : 46 : 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
712 : 46 : write(ab_out,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BV,' and Shear G=',GV
713 : 46 : Eaverage=9.d0*BV*GV/(3*BV+GV)
714 : 46 : Nuaverage=0.5*(1.d0-(3.d0*GV)/(3.d0*BV+GV) )
715 : 46 : Laverage=(3.d0*BV-2.d0*GV)/3.d0
716 : 46 : write(ab_out,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
717 : 92 : & ' and Poisson ratio Nu=',Nuaverage
718 : 46 : Vp=dsqrt(1.d9*(BV+4.d0*GV/3.d0)/rho)
719 : 46 : Vs=dsqrt(1.d9*GV/rho)
720 : 46 : Vphi=dsqrt(1.d9*BV/rho)
721 : 46 : write(ab_out,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
722 : 46 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
723 : 46 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Model%ucvol)**(1./3.)*Ha_K
724 : 46 : write(ab_out,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
725 : :
726 : : ! Reuss notation
727 : 46 : write(ab_out,'(a)')' '
728 : 46 : write(ab_out,'(a)')' ========================= Reuss average (constant stress) ==================='
729 : 46 : BR=1.d0/(Sij(1,1)+Sij(2,2)+Sij(3,3)+2.d0*(Sij(1,2)+Sij(1,3)+Sij(2,3)))
730 : 46 : 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)))
731 : 46 : write(ab_out,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BR,' and Shear G=',GR
732 : 46 : Eaverage=9.d0*BR*GR/(3*BR+GR)
733 : 46 : Nuaverage=0.5*(1.d0-(3.d0*GR)/(3.d0*BR+GR) )
734 : 46 : Laverage=(3.d0*BR-2.d0*GR)/3.d0
735 : 46 : write(ab_out,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
736 : 92 : & ' and Poisson ratio Nu=',Nuaverage
737 : 46 : Vp=dsqrt(1.d9*(BR+4.d0*GR/3.d0)/rho)
738 : 46 : Vs=dsqrt(1.d9*GR/rho)
739 : 46 : Vphi=dsqrt(1.d9*BR/rho)
740 : 46 : write(ab_out,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
741 : 46 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
742 : 46 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Model%ucvol)**(1./3.)*Ha_K
743 : 46 : write(ab_out,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
744 : :
745 : : ! Voigt-Reuss-Hill notation
746 : 46 : write(ab_out,'(a)')' '
747 : 46 : write(ab_out,'(a)')' ============================== Hill average ================================='
748 : 46 : BH=(BR+BV)/2.d0
749 : 46 : GH=(GR+GV)/2.d0
750 : 46 : write(ab_out,'(2(a,f9.3))')' ISOTHERMAL modulus [in GPa]: Bulk Kt=',BH,' and Shear G=',GH
751 : 46 : Eaverage=9.d0*BH*GH/(3*BH+GH)
752 : 46 : Nuaverage=0.5*(1.d0-(3.d0*GH)/(3.d0*BH+GH) )
753 : 46 : Laverage=(3.d0*BH-2.d0*GH)/3.d0
754 : 46 : write(ab_out,'(3(a,f9.3))')' Average of Young modulus E [in GPa]=',Eaverage,' Lame modulus Lambda [in GPa]=',Laverage,&
755 : 92 : & ' and Poisson ratio Nu=',Nuaverage
756 : 46 : Vp=dsqrt(1.d9*(BH+4.d0*GH/3.d0)/rho)
757 : 46 : Vs=dsqrt(1.d9*GH/rho)
758 : 46 : Vphi=dsqrt(1.d9*BH/rho)
759 : 46 : write(ab_out,'(3(a,f9.3,1x))')' Velocities [in m.s-1]: compressional Vp=',Vp,' shear Vs=',Vs,' and bulk Vphi=',Vphi
760 : 46 : V_Debye=(1./3.*(1./Vp**3+2./Vs**3))**(-1./3.)
761 : 46 : T_Debye= (V_Debye/1.d3/Bohr_Ang/1.d-13*Time_Sec)*(six*pi**2*Invar%natom_unitcell/Model%ucvol)**(1./3.)*Ha_K
762 : 46 : write(ab_out,'(2(a,f9.3,1x))')' Debye velocity [in m.s-1]=',V_Debye,' and temperature [in K]=',T_Debye
763 : :
764 : : ! Compute the elastic anisotropy
765 : 46 : A_U=5.d0*GV/GR+BV/BR -6.d0
766 : 46 : A_B=(BV-BR)/(BV+BR)
767 : 46 : A_G=(GV-GR)/(GV+GR)
768 : 46 : write(ab_out,'(a)')' '
769 : 46 : write(ab_out,'(a)')' ========================= Elastic anisotropy ================================='
770 : 46 : write(ab_out,'(a,f9.3)')' Elastic anisotropy index : A_U= 5*G_V/G_R + K_V/K_R - 6 =',A_U
771 : 46 : write(ab_out,'(a,f9.3)')' Bulk anisotropy ratio : A_B= (B_V-B_R)/(B_V+B_R) =',A_B
772 : 46 : write(ab_out,'(a,f9.3)')' Shear anisotropy ratio : A_G= (G_V-G_R)/(G_V+G_R) =',A_G
773 : :
774 : :
775 : : ! Store the RVH value of the elastic moduli (will be useful for the Gruneisen)
776 : 46 : Model%BulkModulus_T=BH
777 : 46 : Model%Shear=GH
778 : 46 : Model%Density=rho
779 : :
780 : 46 : ABI_FREE(cijkl)
781 : 46 : ABI_FREE(Cij)
782 : 46 : ABI_FREE(Sij)
783 : 46 : end subroutine tdep_model_calc_elastic
784 : : !!***
785 : :
786 : : !====================================================================================================
787 : :
788 : : !!****f* m_tdep_model/tdep_model_calc_thermo
789 : : !! NAME
790 : : !! tdep_model_calc_thermo
791 : : !!
792 : : !! FUNCTION
793 : : !! Compute the thermodynamic quantities (heat capacity, entropy, internal
794 : : !! and free energy, mean square displacement, ...) from the phonon DOS,
795 : : !! both directly at Invar%temperature and, under the harmonic
796 : : !! approximation, as a function of temperature. Write the results to the
797 : : !! _thermo.dat file.
798 : : !!
799 : : !! INPUTS
800 : : !! Invar = Input object containing the input variables.
801 : : !! MPIdata = Info on MPI parallelism.
802 : : !! PHdos = Phonon density of states.
803 : : !!
804 : : !! OUTPUT
805 : : !!
806 : : !! SIDE EFFECTS
807 : : !! Model = TDEP Model object. On entry, Model%U0 and Model%ucvol must
808 : : !! already be set. On exit, Model%HeatCapa_V is computed and
809 : : !! stored.
810 : : !!
811 : : !! SOURCE
812 : :
813 : 46 : subroutine tdep_model_calc_thermo(Model,Invar,MPIdata,PHdos)
814 : :
815 : : class(tdep_Model_type), intent(inout) :: Model
816 : : type(atdep_dataset_type),intent(in) :: Invar
817 : : type(MPI_enreg_type), intent(in) :: MPIdata
818 : : type(phdos_t),intent(in) :: PHdos
819 : :
820 : : integer :: iomega,itemp,iatom,itypat
821 : : double precision :: k_B,wovert,heatcapa,entropy,internalE,freeE,expm2x,ln2shx,cothx,xx
822 : : double precision :: Ftot,domega,MSD,Omega_m2,mass_amu,vdos
823 : :
824 : 46 : write(ab_out,'(a)')' '
825 : 46 : write(ab_out,'(a)') ' #############################################################################'
826 : 46 : write(ab_out,'(a)') ' ################# Thermodynamic quantities: Free energy,...##################'
827 : 46 : write(ab_out,'(a)') ' #############################################################################'
828 : :
829 : : ! The heat capacity, entropy, internal and free energies (direct calculation)
830 : : ! ===========================================================================
831 : 46 : mass_amu=zero
832 [ + + ]: 156 : do iatom=1,Invar%natom_unitcell
833 : 110 : itypat=Invar%typat_unitcell(iatom)
834 : 156 : mass_amu=mass_amu+Invar%amu(itypat)
835 : : end do
836 : 46 : mass_amu=mass_amu*amu_emass/real(Invar%natom_unitcell)
837 : :
838 : : ! k_B is Boltzmann constant in eV/K
839 : 46 : k_B=kb_HaK*Ha_eV
840 : 46 : domega=(Invar%dosdeltae*Ha_meV)
841 : 46 : wovert=1.d0/(2*Invar%temperature*k_B)
842 : 46 : heatcapa=0.d0 ; entropy=0.d0 ; internalE=0.d0 ; freeE=0.d0 ; MSD=0.d0 ; Omega_m2=0.d0 ; vdos=0.d0
843 [ + + ]: 90332 : do iomega=1,PHdos%nomega
844 : 90286 : xx=PHdos%omega(iomega)*Ha_eV
845 [ + + ]: 90286 : if (xx.lt.tol8) cycle
846 : 89945 : expm2x=exp(-2.d0*wovert*xx)
847 : 89945 : ln2shx=wovert*xx+log(1.d0-expm2x)
848 : 89945 : cothx=(1.d0+expm2x)/(1.d0-expm2x)
849 : 89945 : heatcapa =heatcapa + (wovert*xx/sinh(wovert*xx))**2*PHdos%phdos(iomega)*domega
850 : 89945 : internalE=internalE + (wovert*xx)*cothx*PHdos%phdos(iomega)*domega
851 : 89945 : entropy =entropy + ((wovert*xx)*cothx-ln2shx)*PHdos%phdos(iomega)*domega
852 : 89945 : freeE =freeE + log(2*sinh(wovert*xx))*PHdos%phdos(iomega)*domega
853 : 89945 : MSD =MSD + (cothx/PHdos%omega(iomega))*PHdos%phdos(iomega)*domega
854 : 89945 : Omega_m2 =Omega_m2 + (1.d0/PHdos%omega(iomega))**2*PHdos%phdos(iomega)*domega
855 : 90332 : vdos =vdos + PHdos%phdos(iomega)*domega
856 : : end do
857 : 46 : heatcapa=heatcapa*3
858 : 46 : entropy=entropy*3
859 : 46 : internalE=internalE*3*k_B*Invar%temperature
860 : 46 : freeE=freeE*3*k_B*Invar%temperature
861 : 46 : MSD=MSD*3.d0/mass_amu/2.d0
862 : 46 : Omega_m2=Omega_m2*3.d0
863 : 46 : Ftot=Model%U0*Ha_eV+freeE
864 : 46 : Model%HeatCapa_V=heatcapa
865 [ - + ]: 46 : if (.not.MPIdata%iam_master) return
866 : : ! End of the calculation --> RETURN
867 : : ! =================================
868 : :
869 : 46 : write(ab_out,'(a)') ' Writing thermodynamical quantities in file thermo.dat'
870 : 46 : open(unit=20,file=trim(Invar%output_prefix)//'_thermo.dat')
871 : 46 : write(20,'(a)')'============= Direct results (without any inter/extrapolation) =================='
872 : 46 : write(20,'(1x,a,f10.3)')'For present temperature (in Kelvin): T= ',Invar%temperature
873 : 46 : write(20,'(1x,a,f12.3)')' The cold contribution (in eV/atom): U_0 =',Model%U0*Ha_eV
874 : 46 : write(20,'(1x,a,f10.3)')' The specific heat (in k_b/atom): C_v=',heatcapa
875 : 46 : write(20,'(1x,a,f10.3)')' The vibrational entropy (in k_b/atom): S_vib =',entropy
876 : 46 : write(20,'(1x,a,f10.3)')' The internal energy (in eV/atom): U_vib =',internalE
877 : 46 : write(20,'(1x,a,f10.3)')' The vibrational contribution (in eV/atom): F_vib = U_vib -T.S_vib =',freeE
878 : 46 : write(20,'(1x,a,f12.3)')' The harmonic free energy (in eV/atom) --> F_tot^HA = U_0 + F_vib =',Ftot
879 : 46 : write(20,'(1x,a)')' Useful quantities for melting :'
880 : 46 : write(20,'(1x,a,f10.3)')' The mean square displacement (in a.u.): sqrt(<u^2>) =',(MSD)**0.5
881 : 46 : write(20,'(1x,a,f10.3)')' The <Omega^(-2)> factor (in THz^(-2)) =',Omega_m2/(Ha_THz)**2
882 : 46 : write(20,'(1x,a,f10.3)')' The Wigner-Seitz radius (in a.u.) : d_at =',(6*Model%ucvol/pi/real(Invar%natom_unitcell))**(1./3.)
883 : 46 : write(20,'(1x,a,f10.3)')' The average mass / proton-electron mass ratio (in a.u.) =', mass_amu/amu_emass
884 : 46 : write(20,'(1x,a,f10.3)')' The Lindemann constant : sqrt(<u^2>)/d_at =',(MSD)**0.5/&
885 : 92 : & (6*Model%ucvol/pi/real(Invar%natom_unitcell))**(1./3.)
886 : 46 : write(20,'(1x,a,f10.3)')' The integral of vDOS =',vdos
887 : 46 : write(20,'(a)')' '
888 : :
889 : : ! The free energy (extrapolation)
890 : : ! ===============================
891 : 46 : write(20,'(a)')'============= Harmonic Approximation (HA) =================='
892 : 46 : write(20,'(1x,a)')' Note that the following results come from an EXTRAPOLATION:'
893 : 46 : write(20,'(1x,a,i5,a)')' 1/ F_vib^HA(T) is computed for each T using vDOS(T=',int(Invar%temperature),')'
894 : 46 : write(20,'(1x,a)')' 2/ F_tot^HA(T) = F_vib^HA(T) + U_0'
895 : : write(20,'(a)')' T F_vib^HA(T) F_tot^HA(T) C_v(T) '&
896 : 46 : & //' S_vib(T) U_vib(T) MSD(T)'
897 [ + + ]: 4646 : do itemp=1,100
898 : 4600 : wovert=1.d0/(2*real(itemp)*100*k_B)
899 : 4600 : freeE=0.d0
900 : 4600 : heatcapa=0.d0
901 : 4600 : entropy=0.d0
902 : 4600 : internalE=0.d0
903 : 4600 : MSD=0.d0
904 [ + + ]: 9033200 : do iomega=1,PHdos%nomega
905 : 9028600 : xx=PHdos%omega(iomega)*Ha_eV
906 [ + + ]: 9028600 : if (xx.lt.tol8) cycle
907 : 8994500 : expm2x=exp(-2.d0*wovert*xx)
908 : 8994500 : ln2shx=wovert*xx+log(1.d0-expm2x)
909 : 8994500 : cothx=(1.d0+expm2x)/(1.d0-expm2x)
910 : 8994500 : heatcapa =heatcapa + (wovert*xx/sinh(wovert*xx))**2*PHdos%phdos(iomega)*domega
911 : 8994500 : internalE=internalE + (wovert*xx)*cothx*PHdos%phdos(iomega)*domega
912 : 8994500 : entropy =entropy + ((wovert*xx)*cothx-ln2shx)*PHdos%phdos(iomega)*domega
913 : 8994500 : freeE =freeE + log(2*sinh(wovert*xx))*PHdos%phdos(iomega)*domega
914 : 9033200 : MSD =MSD + (cothx/PHdos%omega(iomega))*PHdos%phdos(iomega)*domega
915 : : end do
916 : 4600 : heatcapa=heatcapa*3
917 : 4600 : entropy=entropy*3
918 : 4600 : internalE=internalE*3*k_B*itemp*100
919 : 4600 : freeE=freeE*3*k_B*itemp*100
920 : 4600 : MSD=MSD*3.d0/mass_amu/2.d0
921 : 4600 : Ftot=Model%U0*Ha_eV+freeE
922 : 4646 : write(20,'(1x,i5,6(1x,f15.3))') itemp*100,freeE,Ftot,heatcapa,entropy,internalE,(MSD)**0.5
923 : : end do
924 : 46 : close(20)
925 : :
926 : : end subroutine tdep_model_calc_thermo
927 : : !!***
928 : :
929 : : !====================================================================================================
930 : :
931 : : !!****f* m_tdep_model/tdep_calc_ftot2
932 : : !! NAME
933 : : !! tdep_calc_ftot2
934 : : !!
935 : : !! FUNCTION
936 : : !! Compute the second-order IFC contribution to the total forces and energy of the model.
937 : : !! The first-order IFC contribution is also computed here.
938 : : !!
939 : : !! INPUTS
940 : : !! Invar = Input object containing the input variables.
941 : : !! Phi2(3*natom,3*natom) = Second-order IFC, in the supercell.
942 : : !! ucart(3,natom,my_nstep) = Cartesian displacements of the atoms at each MD step.
943 : : !!
944 : : !! OUTPUT
945 : : !!
946 : : !! SIDE EFFECTS
947 : : !! Model = TDEP Model object. On entry, Model%Phi1 must already be set.
948 : : !! On exit, Model%Phi1Ui and Model%Phi2UiUj are computed, and
949 : : !! Model%Forces is overwritten with the first- and second-order
950 : : !! contributions: F = -Phi1 - Phi2.u
951 : : !!
952 : : !! SOURCE
953 : :
954 : 51 : subroutine tdep_calc_ftot2(Model,Invar,Phi2,ucart)
955 : :
956 : : class(tdep_Model_type),intent(inout) :: Model
957 : : type(atdep_dataset_type),intent(in) :: Invar
958 : : double precision, intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
959 : : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
960 : :
961 : : integer :: jj,istep,jatom
962 : 51 : double precision, allocatable :: ucart_blas(:)
963 : 51 : double precision, allocatable :: ftot2(:)
964 : :
965 : : ! Compute Forces of the model TDEP
966 [ + - - + : 24201 : ABI_CALLOC(ucart_blas,(3*Invar%natom))
+ - + - +
- - + + -
+ + ]
967 [ - + + - : 24099 : ABI_CALLOC(ftot2 ,(3*Invar%natom))
- + + - +
+ ]
968 [ + + ]: 1143 : do istep=1,Invar%my_nstep
969 [ + + ]: 484548 : ucart_blas(:)=0.d0
970 [ + + ]: 484548 : ftot2 (:)=0.d0
971 [ + + ]: 162244 : do jatom=1,Invar%natom
972 [ + + ]: 645700 : do jj=1,3
973 : 644608 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
974 : : end do
975 : : end do
976 [ + + ]: 484548 : Model%Phi1Ui(istep)=sum(Model%Phi1(:)*ucart_blas(:))
977 : 1092 : call DGEMM('N','N',3*Invar%natom,1,3*Invar%natom,1.d0,Phi2,3*Invar%natom,ucart_blas,3*Invar%natom,0.d0,ftot2(:),3*Invar%natom)
978 : 1092 : call DGEMM('T','N',1,1,3*Invar%natom,1./2.d0,ftot2,3*Invar%natom,ucart_blas,3*Invar%natom,0.d0,Model%Phi2UiUj(istep),3*Invar%natom)
979 [ + + ]: 484599 : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=-Model%Phi1(:)-ftot2(:)
980 : : end do !istep
981 [ - + ]: 51 : ABI_FREE(ucart_blas)
982 [ - + ]: 51 : ABI_FREE(ftot2)
983 : :
984 : 51 : end subroutine tdep_calc_ftot2
985 : : !!***
986 : :
987 : : !====================================================================================================
988 : :
989 : : !!****f* m_tdep_model/tdep_calc_ftot3
990 : : !! NAME
991 : : !! tdep_calc_ftot3
992 : : !!
993 : : !! FUNCTION
994 : : !! Compute the third-order IFC contribution to the total forces and energy of the model.
995 : : !!
996 : : !! INPUTS
997 : : !! Invar = Input object containing the input variables.
998 : : !! Shell3at = Shell object for third-order interactions.
999 : : !! ucart(3,natom,my_nstep) = Cartesian displacements of the atoms at each MD step.
1000 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1001 : : !!
1002 : : !! OUTPUT
1003 : : !!
1004 : : !! SIDE EFFECTS
1005 : : !! Model = TDEP Model object. On entry, Model%Phi3 must already be set,
1006 : : !! and Model%Forces must contain the lower-order contributions
1007 : : !! (see tdep_calc_ftot2). On exit, Model%Phi3UiUjUk is computed,
1008 : : !! and the third-order contribution, -1/2 Phi3.u.u, is added to
1009 : : !! Model%Forces.
1010 : : !!
1011 : : !! SOURCE
1012 : :
1013 : 16 : subroutine tdep_calc_ftot3(Model,Invar,Shell3at,ucart,Sym)
1014 : :
1015 : : class(tdep_Model_type),intent(inout) :: Model
1016 : : type(atdep_dataset_type),intent(in) :: Invar
1017 : : type(Shell_type),intent(in) :: Shell3at
1018 : : type(Symmetries_type),intent(in) :: Sym
1019 : : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
1020 : :
1021 : : integer :: iatom,jatom,katom,isym,itrans,ishell,iatshell
1022 : : integer :: ii,jj,kk,istep
1023 : 16 : double precision, allocatable :: Phi3_333(:,:,:)
1024 : 16 : double precision, allocatable :: ucart_blas(:)
1025 : 16 : double precision, allocatable :: ftot3(:,:)
1026 : :
1027 [ + - + + : 656 : ABI_MALLOC(Phi3_333,(3,3,3)) ; Phi3_333(:,:,:)=0.d0
+ + + + ]
1028 [ + - - + : 176560 : ABI_MALLOC(ftot3,(3*Invar%natom,Invar%my_nstep)); ftot3(:,:)=0.d0
+ - - + +
- + - + -
- + - + +
+ + + ]
1029 [ + + ]: 2952 : do iatom=1,Invar%natom
1030 [ + + ]: 12792 : do ishell=1,Shell3at%nshell
1031 : : ! Build the 3x3x3 IFC of an atom in this shell
1032 [ + + ]: 9840 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1033 [ + + ]: 108992 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
1034 : 97496 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1035 : 97496 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1036 : 97496 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1037 : 97496 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1038 : 97496 : call tdep_build_phi3_333(isym,Model%Phi3(:,:,:,ishell),Phi3_333,Sym,itrans)
1039 : : ! Calculation of the force components (third order)
1040 [ + + ]: 2057256 : do istep=1,Invar%my_nstep
1041 [ + + ]: 7897176 : do ii=1,3
1042 [ + + ]: 25348960 : do jj=1,3
1043 [ + + ]: 76046880 : do kk=1,3
1044 : : ftot3(3*(iatom-1)+ii,istep)=ftot3(3*(iatom-1)+ii,istep)+&
1045 : 70197120 : & Phi3_333(ii,jj,kk)*ucart(jj,jatom,istep)*ucart(kk,katom,istep)
1046 : : end do !kk
1047 : : end do !jj
1048 : : end do !ii
1049 : : end do !istep
1050 : : end do !iatshell
1051 : : end do !ishell
1052 : : end do !iatom
1053 [ - + ]: 16 : ABI_FREE(Phi3_333)
1054 [ + + + + ]: 176496 : ftot3(:,:)=ftot3(:,:)/2.d0
1055 : :
1056 [ - + + - : 8872 : ABI_MALLOC(ucart_blas,(3*Invar%natom)); ucart_blas(:)=0.d0
+ - + - -
+ + - +
+ ]
1057 [ + + ]: 336 : do istep=1,Invar%my_nstep
1058 [ + + ]: 176480 : ucart_blas(:)=0.d0
1059 [ + + ]: 59040 : do jatom=1,Invar%natom
1060 [ + + ]: 235200 : do jj=1,3
1061 : 234880 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
1062 : : end do
1063 : : end do
1064 : : call DGEMM('T','N',1,1,3*Invar%natom,1./3.d0,ftot3(:,istep),3*Invar%natom,ucart_blas,&
1065 : 320 : & 3*Invar%natom,0.d0,Model%Phi3UiUjUk(istep),3*Invar%natom)
1066 : : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=&
1067 [ + + ]: 176496 : & Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)-ftot3(:,istep)
1068 : : end do
1069 [ - + ]: 16 : ABI_FREE(ucart_blas)
1070 [ - + ]: 16 : ABI_FREE(ftot3)
1071 : :
1072 : 16 : end subroutine tdep_calc_ftot3
1073 : : !!***
1074 : :
1075 : : !====================================================================================================
1076 : :
1077 : : !!****f* m_tdep_model/tdep_calc_ftot4
1078 : : !! NAME
1079 : : !! tdep_calc_ftot4
1080 : : !!
1081 : : !! FUNCTION
1082 : : !! Compute the fourth-order IFC contribution to the total forces and energy of the model.
1083 : : !!
1084 : : !! INPUTS
1085 : : !! Invar = Input object containing the input variables.
1086 : : !! Shell4at = Shell object for fourth-order interactions.
1087 : : !! ucart(3,natom,my_nstep) = Cartesian displacements of the atoms at each MD step.
1088 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1089 : : !!
1090 : : !! OUTPUT
1091 : : !!
1092 : : !! SIDE EFFECTS
1093 : : !! Model = TDEP Model object. On entry, Model%Phi4 must already be set,
1094 : : !! and Model%Forces must contain the lower-order contributions
1095 : : !! (see tdep_calc_ftot2 and tdep_calc_ftot3). On exit,
1096 : : !! Model%Phi4UiUjUkUl is computed, and the fourth-order
1097 : : !! contribution, -1/6 Phi4.u.u.u, is added to Model%Forces.
1098 : : !!
1099 : : !! SOURCE
1100 : :
1101 : 9 : subroutine tdep_calc_ftot4(Model,Invar,Shell4at,ucart,Sym)
1102 : :
1103 : : class(tdep_Model_type),intent(inout) :: Model
1104 : : type(atdep_dataset_type),intent(in) :: Invar
1105 : : type(Shell_type),intent(in) :: Shell4at
1106 : : type(Symmetries_type),intent(in) :: Sym
1107 : : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
1108 : :
1109 : : integer :: iatom,jatom,katom,latom,isym,itrans,ishell,iatshell
1110 : : integer :: ii,jj,kk,ll,istep
1111 : 9 : double precision, allocatable :: Phi4_3333(:,:,:,:)
1112 : 9 : double precision, allocatable :: ucart_blas(:)
1113 : 9 : double precision, allocatable :: ftot4(:,:)
1114 : :
1115 [ + - + + : 1098 : ABI_MALLOC(Phi4_3333,(3,3,3,3)) ; Phi4_3333(:,:,:,:)=0.d0
+ + + + +
+ ]
1116 [ + - - + : 101745 : ABI_MALLOC(ftot4,(3*Invar%natom,Invar%my_nstep)); ftot4(:,:)=0.d0
+ - - + +
- + - + -
- + - + +
+ + + ]
1117 [ + + ]: 1701 : do iatom=1,Invar%natom
1118 [ + + ]: 10181 : do ishell=1,Shell4at%nshell
1119 : : ! Build the 3x3x3x3 IFC of an atom in this shell
1120 [ + + ]: 8480 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1121 [ + + ]: 259144 : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1122 : 249996 : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1123 : 249996 : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1124 : 249996 : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
1125 : 249996 : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1126 : 249996 : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1127 : 249996 : call tdep_build_phi4_3333(isym,Model%Phi4(:,:,:,:,ishell),Phi4_3333,Sym,itrans)
1128 : : ! Calculation of the force components (third order)
1129 [ + + ]: 5258396 : do istep=1,Invar%my_nstep
1130 [ + + ]: 20249676 : do ii=1,3
1131 [ + + ]: 64998960 : do jj=1,3
1132 [ + + ]: 194996880 : do kk=1,3
1133 [ + + ]: 584990640 : do ll=1,3
1134 : : ftot4(3*(iatom-1)+ii,istep)=ftot4(3*(iatom-1)+ii,istep)+&
1135 : 539991360 : & Phi4_3333(ii,jj,kk,ll)*ucart(jj,jatom,istep)*ucart(kk,katom,istep)*ucart(ll,latom,istep)
1136 : : end do !ll
1137 : : end do !kk
1138 : : end do !jj
1139 : : end do !ii
1140 : : end do !istep
1141 : : end do !iatshell
1142 : : end do !ishell
1143 : : end do !iatom
1144 [ - + ]: 9 : ABI_FREE(Phi4_3333)
1145 [ + + + + ]: 101709 : ftot4(:,:)=ftot4(:,:)/6.d0
1146 : :
1147 [ - + + - : 5112 : ABI_MALLOC(ucart_blas ,(3*Invar%natom)) ; ucart_blas (:)=0.d0
+ - + - -
+ + - +
+ ]
1148 [ + + ]: 189 : do istep=1,Invar%my_nstep
1149 [ + + ]: 101700 : ucart_blas(:)=0.d0
1150 [ + + ]: 34020 : do jatom=1,Invar%natom
1151 [ + + ]: 135540 : do jj=1,3
1152 : 135360 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
1153 : : end do
1154 : : end do
1155 : : call DGEMM('T','N',1,1,3*Invar%natom,1./4.d0,ftot4(:,istep),3*Invar%natom,ucart_blas,&
1156 : 180 : & 3*Invar%natom,0.d0,Model%Phi4UiUjUkUl(istep),3*Invar%natom)
1157 : : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=&
1158 [ + + ]: 101709 : & Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)-ftot4(:,istep)
1159 : : end do
1160 [ - + ]: 9 : ABI_FREE(ucart_blas)
1161 [ - + ]: 9 : ABI_FREE(ftot4)
1162 : :
1163 : 9 : end subroutine tdep_calc_ftot4
1164 : : !!***
1165 : :
1166 : : !====================================================================================================
1167 : :
1168 : : !!****f* m_tdep_model/tdep_model_write_phi1
1169 : : !! NAME
1170 : : !! tdep_model_write_phi1
1171 : : !!
1172 : : !! FUNCTION
1173 : : !! Write the first-order IFC of the model in the output file
1174 : : !! (see tdep_write_phi1 in m_tdep_phi2).
1175 : : !!
1176 : : !! INPUTS
1177 : : !! Model = TDEP Model object, containing the first-order IFC (Model%Phi1).
1178 : : !!
1179 : : !! OUTPUT
1180 : : !! Only writing.
1181 : : !!
1182 : : !! SOURCE
1183 : :
1184 : 46 : subroutine tdep_model_write_phi1(Model)
1185 : :
1186 : : class(tdep_Model_type),intent(in) :: Model
1187 : :
1188 : 46 : call tdep_write_phi1(Model%Phi1, Model%natom, Model%natom_unitcell)
1189 : :
1190 : 46 : end subroutine tdep_model_write_phi1
1191 : : !!***
1192 : :
1193 : : !====================================================================================================
1194 : :
1195 : : !!****f* m_tdep_model/tdep_model_write_phi2
1196 : : !! NAME
1197 : : !! tdep_model_write_phi2
1198 : : !!
1199 : : !! FUNCTION
1200 : : !! Write the (short-range) second-order IFC of the model in the output file,
1201 : : !! for each shell, and in the _Phi2*.dat files
1202 : : !! (see tdep_write_phi2 in m_tdep_phi2).
1203 : : !!
1204 : : !! INPUTS
1205 : : !! Model = TDEP Model object, containing the second-order IFC (Model%Phi2%SR).
1206 : : !! Invar = Input object containing the input variables.
1207 : : !! MPIdata = Info on MPI parallelism.
1208 : : !! SC = Supercell object describing the ideal structure.
1209 : : !! Shell2at = Shell object for second-order interactions.
1210 : : !!
1211 : : !! OUTPUT
1212 : : !! Only writing.
1213 : : !!
1214 : : !! SOURCE
1215 : :
1216 : 46 : subroutine tdep_model_write_phi2(Model,Invar,MPIdata,SC,Shell2at)
1217 : :
1218 : : class(tdep_Model_type),intent(in) :: Model
1219 : : type(atdep_dataset_type),intent(in) :: Invar
1220 : : type(MPI_enreg_type), intent(in) :: MPIdata
1221 : : type(tdep_Supercell_type),intent(in) :: SC
1222 : : type(Shell_type),intent(in) :: Shell2at
1223 : :
1224 : 46 : call tdep_write_phi2(Invar,MPIdata,SC,Model%Phi2%SR,Shell2at)
1225 : :
1226 : 46 : end subroutine tdep_model_write_phi2
1227 : : !!***
1228 : :
1229 : : !====================================================================================================
1230 : :
1231 : : !!****f* m_tdep_model/tdep_model_write_phi4
1232 : : !! NAME
1233 : : !! tdep_model_write_phi4
1234 : : !!
1235 : : !! FUNCTION
1236 : : !! Write the fourth-order IFC of the model in the output file, for each shell
1237 : : !! (see tdep_write_phi4 in m_tdep_phi4).
1238 : : !!
1239 : : !! INPUTS
1240 : : !! Model = TDEP Model object, containing the fourth-order IFC (Model%Phi4).
1241 : : !! SC = Supercell object describing the ideal structure.
1242 : : !! Shell4at = Shell object for fourth-order interactions.
1243 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1244 : : !!
1245 : : !! OUTPUT
1246 : : !! Only writing.
1247 : : !!
1248 : : !! SOURCE
1249 : :
1250 : 9 : subroutine tdep_model_write_phi4(Model,SC,Shell4at,Sym)
1251 : :
1252 : : class(tdep_Model_type),intent(in) :: Model
1253 : : type(tdep_Supercell_type),intent(in) :: SC
1254 : : type(Shell_type),intent(in) :: Shell4at
1255 : : type(Symmetries_type),intent(in) :: Sym
1256 : :
1257 : 9 : call tdep_write_phi4(Model%Phi4,SC,Shell4at,Sym)
1258 : :
1259 : 9 : end subroutine tdep_model_write_phi4
1260 : : !!***
1261 : :
1262 : : !====================================================================================================
1263 : :
1264 : : !!****f* m_tdep_model/tdep_model_write_phi3
1265 : : !! NAME
1266 : : !! tdep_model_write_phi3
1267 : : !!
1268 : : !! FUNCTION
1269 : : !! Write the third-order IFC of the model in the output file, for each shell
1270 : : !! (see tdep_write_phi3 in m_tdep_phi3).
1271 : : !!
1272 : : !! INPUTS
1273 : : !! Model = TDEP Model object, containing the third-order IFC (Model%Phi3).
1274 : : !! SC = Supercell object describing the ideal structure.
1275 : : !! Shell3at = Shell object for third-order interactions.
1276 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1277 : : !!
1278 : : !! OUTPUT
1279 : : !! Only writing.
1280 : : !!
1281 : : !! SOURCE
1282 : :
1283 : 16 : subroutine tdep_model_write_phi3(Model,SC,Shell3at,Sym)
1284 : :
1285 : : class(tdep_Model_type),intent(inout) :: Model
1286 : : type(tdep_Supercell_type),intent(in) :: SC
1287 : : type(Shell_type),intent(in) :: Shell3at
1288 : : type(Symmetries_type),intent(in) :: Sym
1289 : :
1290 : 16 : call tdep_write_phi3(Model%Phi3,SC,Shell3at,Sym)
1291 : :
1292 : 16 : end subroutine tdep_model_write_phi3
1293 : : !!***
1294 : :
1295 : : !====================================================================================================
1296 : :
1297 : : !!****f* m_tdep_model/tdep_check_constraints
1298 : : !! NAME
1299 : : !! tdep_check_constraints
1300 : : !!
1301 : : !! FUNCTION
1302 : : !! Check that the IFC of the model fulfill the acoustic sum rules and the
1303 : : !! invariance under an arbitrary rotation of the system, for the first,
1304 : : !! second and third orders (the checks at the fourth order are currently
1305 : : !! disabled). A warning is issued for each constraint that is not fulfilled.
1306 : : !!
1307 : : !! INPUTS
1308 : : !! Model = TDEP Model object, containing the IFC at all orders.
1309 : : !! Invar = Input object containing the input variables.
1310 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
1311 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1312 : : !! Shell3at = Shell object for third-order interactions.
1313 : : !! Shell4at = Shell object for fourth-order interactions.
1314 : : !!
1315 : : !! OUTPUT
1316 : : !! Only warnings are printed.
1317 : : !!
1318 : : !! SOURCE
1319 : :
1320 : 46 : subroutine tdep_check_constraints(Model,Invar,SC,Sym,Shell3at,Shell4at)
1321 : :
1322 : : class(tdep_Model_type),intent(in) :: Model
1323 : : type(tdep_Supercell_type),intent(in) :: SC
1324 : : type(atdep_dataset_type),intent(in) :: Invar
1325 : : type(Symmetries_type),intent(in) :: Sym
1326 : : type(Shell_type), intent(in) :: Shell3at
1327 : : type(Shell_type), intent(in) :: Shell4at
1328 : :
1329 : : integer :: natom,natom_unitcell
1330 : : integer :: ii,jj,kk,ll,iatom,jatom,katom,latom,isym,itrans
1331 : : integer :: alpha,beta,gama,lambda
1332 : : integer :: ishell,iatshell
1333 : : double precision :: norm1
1334 : : double precision :: Kroenecker(3,3),Phi3_333(3,3,3),Phi4_3333(3,3,3,3)
1335 : 46 : double precision, allocatable :: asr3(:,:,:,:,:),rot3(:,:,:,:,:)
1336 : 46 : double precision, allocatable :: asr4(:,:,:,:,:,:)
1337 : : integer :: ierr
1338 : : logical :: order2,order3,order4
1339 : :
1340 : 46 : ierr = 0;
1341 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Compute the acoustic sum rule and the !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343 : : !!!!!!!!!!!!!!!!!!!!!!!!!!! invariance under an arbitrary rotation of the system !!!!!!!!!!!!!!!!!!!!!!!
1344 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1345 : 46 : natom = Model%natom
1346 : 46 : natom_unitcell = Model%natom_unitcell
1347 : :
1348 : 46 : Kroenecker(:,:)=zero
1349 : 46 : Kroenecker(1,1)=1 ; Kroenecker(2,2)=1 ; Kroenecker(3,3)=1
1350 : :
1351 : 46 : order2 = .false.
1352 : 46 : order3 = .false.
1353 : 46 : order4 = .false.
1354 [ + - ]: 46 : if (Model%order.ge.2) order2=.true.
1355 [ + + ]: 46 : if (Model%order.ge.3) order3=.true.
1356 : : !FB4th if (Model%order.ge.4) order4=.true.
1357 : : if (Model%order.ge.4) order4=.false.
1358 : :
1359 [ + - ]: 46 : if (order2) then
1360 : : ! FIRST ORDER
1361 : : ! Acoustic sum rule (first order)'
1362 [ + + ]: 184 : do alpha=1,3
1363 : 138 : norm1=zero
1364 [ + + ]: 468 : do iatom=1,natom_unitcell
1365 : 468 : norm1=norm1+Model%Phi1((iatom-1)*3+alpha)
1366 : : enddo
1367 [ - + ]: 184 : if (abs(norm1).gt.tol8) then
1368 : 0 : write(std_out,'(a,1(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ASR1) : alpha,Sum =',alpha,norm1
1369 [ # # ]: 0 : if (abs(norm1).gt.tol6) then
1370 : 0 : ABI_WARNING('The acoustic sum rule is not fulfilled (order 1)')
1371 : : end if
1372 : : end if
1373 : : enddo
1374 : : ! Invariance under an arbitrary rotation (first order)'
1375 [ + + ]: 184 : do alpha=1,3
1376 [ + + ]: 598 : do beta=1,3
1377 : 414 : norm1=zero
1378 [ + + ]: 1404 : do iatom=1,natom_unitcell
1379 : : norm1=norm1+Model%Phi1(3*(iatom-1)+alpha)*SC%dij(1,iatom,beta +1)-&
1380 : 1404 : & Model%Phi1(3*(iatom-1)+beta )*SC%dij(1,iatom,alpha+1)
1381 : : end do !iatom
1382 [ - + ]: 552 : if (abs(norm1).gt.tol8) then
1383 : 0 : write(std_out,'(a,2(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ROT1) : alpha,beta,Sum=',alpha,beta,norm1
1384 [ # # ]: 0 : if (abs(norm1).gt.tol6) then
1385 : 0 : ABI_WARNING('The invariance under arbitrary rotation is not fulfilled (order 1)')
1386 : : end if
1387 : : end if
1388 : : end do !beta
1389 : : end do !alpha
1390 : :
1391 : : ! SECOND ORDER
1392 : : ! Acoustic sum rule (second order)
1393 [ + + ]: 7324 : do iatom=1,natom
1394 [ + + ]: 29158 : do alpha=1,3
1395 [ + + ]: 94614 : do beta=1,3
1396 : 65502 : norm1=zero
1397 [ + + ]: 14027778 : do jatom=1,natom
1398 : 14027778 : norm1=norm1+Model%Phi2%SR((iatom-1)*3+alpha,3*(jatom-1)+beta)
1399 : : enddo
1400 [ - + ]: 87336 : if (abs(norm1).gt.tol8) then
1401 : 0 : write(std_out,'(a,3(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ASR2) : iatom,alpha,beta,Sum =',iatom,alpha,beta,norm1
1402 [ # # ]: 0 : if (abs(norm1).gt.tol6) then
1403 : 0 : ABI_WARNING('The acoustic sum rule is not fulfilled (order 2)')
1404 : : end if
1405 : : end if
1406 : : enddo
1407 : : enddo
1408 : : enddo
1409 : : ! Invariance under an arbitrary rotation (first and second order)'
1410 [ + + ]: 7324 : do iatom=1,natom
1411 [ + + ]: 29158 : do alpha=1,3
1412 [ + + ]: 94614 : do beta=1,3
1413 [ + + ]: 283842 : do gama=1,3
1414 : 196506 : norm1=zero
1415 [ + + ]: 42083334 : do jatom=1,natom
1416 : : norm1=norm1+Model%Phi2%SR(3*(iatom-1)+alpha,3*(jatom-1)+beta)*SC%dij(iatom,jatom,gama+1)-&
1417 : 42083334 : & Model%Phi2%SR(3*(iatom-1)+alpha,3*(jatom-1)+gama)*SC%dij(iatom,jatom,beta+1)
1418 : : end do !jatom
1419 : : norm1=norm1+Model%Phi1(3*(iatom-1)+beta)*Kroenecker(alpha,gama)&
1420 : 196506 : & -Model%Phi1(3*(iatom-1)+gama)*Kroenecker(alpha,beta)
1421 [ + + ]: 262008 : if (abs(norm1).gt.tol8) then
1422 : 9184 : write(std_out,'(a,4(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ROT2) : iatom,alpha,beta,gama,Sum =',&
1423 : 18368 : & iatom,alpha,beta,gama,norm1
1424 [ - + ]: 9184 : if (abs(norm1).gt.tol4) then
1425 : 0 : ABI_WARNING('The invariance under arbitrary rotation is not fulfilled (order 2)')
1426 : : end if
1427 : : end if
1428 : : end do !gama
1429 : : end do !beta
1430 : : end do !alpha
1431 : : end do !iatom
1432 : : end if !order=1,2
1433 : :
1434 [ + + ]: 46 : if (order3) then
1435 : : ! THIRD ORDER
1436 [ - + + - : 239832 : ABI_MALLOC(rot3,(natom,3,3,3,3)) ; rot3(:,:,:,:,:)=0.d0
+ - + - +
- + - + -
- + - + +
+ + + + +
+ + + + ]
1437 [ - + + - : 238552 : ABI_MALLOC(asr3,(2,natom,3,3,3)) ; asr3(:,:,:,:,:)=0.d0
+ - + - +
- + - + -
- + + - +
+ + + + +
+ + + + ]
1438 [ + + ]: 2952 : do iatom=1,natom
1439 [ + + + + : 46713176 : rot3(:,:,:,:,:)=0.d0
+ + + + +
+ ]
1440 [ + + + + : 46475360 : asr3(:,:,:,:,:)=0.d0
+ + + + +
+ ]
1441 [ + + ]: 575256 : do jatom=1,natom
1442 [ + + ]: 572320 : if (SC%dij(iatom,jatom,1).gt.Invar%rcut3) cycle
1443 : : ! Compute the rotational invariance (third order)
1444 [ + + ]: 74712 : do alpha=1,3
1445 [ + + ]: 787648 : do beta=1,3
1446 [ + + ]: 699816 : do gama=1,3
1447 [ + + ]: 2099448 : do lambda=1,3
1448 : : rot3(jatom,alpha,beta,gama,lambda)=rot3(jatom,alpha,beta,gama,lambda)+&
1449 : : & Model%Phi2%SR(3*(iatom-1)+gama ,3*(jatom-1)+beta )*Kroenecker(alpha,lambda)+&
1450 : : & Model%Phi2%SR(3*(iatom-1)+alpha ,3*(jatom-1)+gama )*Kroenecker(beta,lambda)-&
1451 : : & Model%Phi2%SR(3*(iatom-1)+lambda,3*(jatom-1)+beta )*Kroenecker(alpha,gama)-&
1452 : 1937952 : & Model%Phi2%SR(3*(iatom-1)+alpha ,3*(jatom-1)+lambda)*Kroenecker(beta,gama)
1453 : : end do !lambda
1454 : : end do !gama
1455 : : end do !beta
1456 : : end do !alpha
1457 : : end do !jatom
1458 [ + + ]: 12776 : do ishell=1,Shell3at%nshell
1459 : : ! Build the 3x3x3 IFC of an atom in this shell
1460 [ + + ]: 9840 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1461 [ + + ]: 108992 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
1462 : 97496 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1463 : 97496 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1464 : 97496 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1465 : 97496 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1466 : 97496 : call tdep_build_phi3_333(isym,Model%Phi3(:,:,:,ishell),Phi3_333,Sym,itrans)
1467 : : ! Compute the first ASR : sum_k Phi3_ijk=0
1468 : : ! --> Phi3_iji+sum_{k.ne.i} Phi3_ijk=0
1469 : : ! --> if i.eq.j Phi3_iii+sum_{k.ne.i} Phi3_iik
1470 [ + + + + : 3899840 : asr3(1,jatom,:,:,:)=asr3(1,jatom,:,:,:)+Phi3_333(:,:,:)
+ + ]
1471 : : ! Compute the second ASR : sum_j Phi3_ijk=0
1472 [ + + + + : 3899840 : asr3(2,katom,:,:,:)=asr3(2,katom,:,:,:)+Phi3_333(:,:,:)
+ + ]
1473 : : ! Compute the rotational invariance (third order)
1474 [ + + ]: 399824 : do alpha=1,3
1475 [ + + ]: 1267448 : do beta=1,3
1476 [ + + ]: 3802344 : do gama=1,3
1477 [ + + ]: 11407032 : do lambda=1,3
1478 : : rot3(jatom,alpha,beta,gama,lambda)=rot3(jatom,alpha,beta,gama,lambda)+&
1479 : : & Phi3_333(alpha,beta,gama )*SC%dij(iatom,katom,lambda+1)-&
1480 : 10529568 : & Phi3_333(alpha,beta,lambda)*SC%dij(iatom,katom,gama +1)
1481 : : end do
1482 : : end do
1483 : : end do
1484 : : end do
1485 : : end do !iatshell
1486 : : end do !ishell
1487 : : ! Check the acoustic sum rules (third order)
1488 [ + + ]: 11744 : do ii=1,3
1489 [ + + ]: 38168 : do jj=1,3
1490 [ + + ]: 114504 : do kk=1,3
1491 : : ! Check the first acoustic sum rule
1492 [ + + ]: 15531912 : do jatom=1,natom
1493 [ - + ]: 15531912 : if (abs(asr3(1,jatom,ii,jj,kk)).gt.tol8) then
1494 : 0 : write(std_out,'(a,1x,5(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ASR3) --->',&
1495 : 0 : & ii,jj,kk,iatom,jatom,asr3(1,jatom,ii,jj,kk)
1496 [ # # ]: 0 : if (abs(asr3(1,jatom,ii,jj,kk)).gt.tol6)&
1497 : 0 : & ABI_WARNING('The acoustic sum rule is not fulfilled (order 3, 3rd dim)')
1498 : : end if
1499 : : end do !jatom
1500 : : ! Check the second acoustic sum rule
1501 [ + + ]: 15558336 : do katom=1,natom
1502 [ - + ]: 15531912 : if (abs(asr3(2,katom,ii,jj,kk)).gt.tol8) then
1503 : 0 : write(std_out,'(a,1x,5(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ASR3) --->',&
1504 : 0 : & ii,jj,kk,iatom,katom,asr3(2,katom,ii,jj,kk)
1505 [ # # ]: 0 : if (abs(asr3(2,katom,ii,jj,kk)).gt.tol6)&
1506 : 0 : & ABI_WARNING('The acoustic sum rule is not fulfilled (order 3, 2nd dim)')
1507 : : end if
1508 : : end do !jatom
1509 : : end do !kk
1510 : : end do !jj
1511 : : end do !ii
1512 : : ! Check the rotational invariance (third order)
1513 [ + + ]: 575272 : do jatom=1,natom
1514 [ + + ]: 2292216 : do alpha=1,3
1515 [ + + ]: 7440160 : do beta=1,3
1516 [ + + ]: 22320480 : do gama=1,3
1517 [ + + ]: 66961440 : do lambda=1,3
1518 [ - + ]: 61810560 : if (abs(rot3(jatom,alpha,beta,gama,lambda)).gt.tol8) then
1519 : 0 : write(std_out,'(a,6(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ROT3) ---> iatom,jatom,alpha,beta,gama,lambda,norm =',&
1520 : 0 : & iatom,jatom,alpha,beta,gama,lambda,rot3(jatom,alpha,beta,gama,lambda)
1521 [ # # ]: 0 : if (abs(rot3(jatom,alpha,beta,gama,lambda)).gt.tol6)&
1522 : 0 : & ABI_WARNING('The invariance under arbitrary rotation is not fulfilled (order 3)')
1523 : : end if
1524 : : end do !lambda
1525 : : end do !gama
1526 : : end do !beta
1527 : : end do !alpha
1528 : : end do !jatom
1529 : : end do !iatom
1530 [ - + ]: 16 : ABI_FREE(asr3)
1531 [ - + ]: 16 : ABI_FREE(rot3)
1532 : : end if !order=3
1533 : :
1534 : : ! FOURTH ORDER
1535 : : if (order4) then
1536 : : !FB ABI_MALLOC(rot3,(natom,3,3,3,3)) ; rot3(:,:,:,:,:)=0.d0
1537 : : ABI_MALLOC(asr4,(natom,natom,3,3,3,3)) ; asr4(:,:,:,:,:,:)=0.d0
1538 : : do iatom=1,natom
1539 : : asr4(:,:,:,:,:,:)=0.d0
1540 : : !FB rot3(:,:,:,:,:)=0.d0
1541 : : !FB do jatom=1,natom
1542 : : !FB if (SC%dij(iatom,jatom,1).gt.rcut3) cycle
1543 : : !FB! Compute the rotational invariance (third order)
1544 : : !FB do alpha=1,3
1545 : : !FB do beta=1,3
1546 : : !FB do gama=1,3
1547 : : !FB do lambda=1,3
1548 : : !FB rot3(jatom,alpha,beta,gama,lambda)=rot3(jatom,alpha,beta,gama,lambda)+&
1549 : : !FB& Model%Phi2%SR(3*(iatom-1)+gama ,3*(jatom-1)+beta )*Kroenecker(alpha,lambda)+&
1550 : : !FB& Model%Phi2%SR(3*(iatom-1)+alpha ,3*(jatom-1)+gama )*Kroenecker(beta,lambda)-&
1551 : : !FB& Model%Phi2%SR(3*(iatom-1)+lambda,3*(jatom-1)+beta )*Kroenecker(alpha,gama)-&
1552 : : !FB& Model%Phi2%SR(3*(iatom-1)+alpha ,3*(jatom-1)+lambda)*Kroenecker(beta,gama)
1553 : : !FB end do !lambda
1554 : : !FB end do !gama
1555 : : !FB end do !beta
1556 : : !FB end do !alpha
1557 : : !FB end do !jatom
1558 : : do ishell=1,Shell4at%nshell
1559 : : ! Build the 3x3x3x3 IFC of an atom in this shell
1560 : : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1561 : : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1562 : : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1563 : : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1564 : : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
1565 : : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1566 : : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1567 : : call tdep_build_phi4_3333(isym,Model%Phi4(:,:,:,:,ishell),Phi4_3333,Sym,itrans)
1568 : : ! Compute the first ASR : sum_l Phi3_ijkl=0
1569 : : asr4(jatom,katom,:,:,:,:)=asr4(jatom,katom,:,:,:,:)+Phi4_3333(:,:,:,:)
1570 : : !FB! Compute the rotational invariance (third order)
1571 : : !FB do alpha=1,3
1572 : : !FB do beta=1,3
1573 : : !FB do gama=1,3
1574 : : !FB do lambda=1,3
1575 : : !FB rot3(jatom,alpha,beta,gama,lambda)=rot3(jatom,alpha,beta,gama,lambda)+&
1576 : : !FB& Phi3_333(alpha,beta,gama )*SC%dij(iatom,katom,lambda+1)-&
1577 : : !FB& Phi3_333(alpha,beta,lambda)*SC%dij(iatom,katom,gama +1)
1578 : : !FB end do
1579 : : !FB end do
1580 : : !FB end do
1581 : : !FB end do
1582 : : end do !iatshell
1583 : : end do !ishell
1584 : : ! Check the acoustic sum rules (fourth order)
1585 : : do ii=1,3
1586 : : do jj=1,3
1587 : : do kk=1,3
1588 : : do ll=1,3
1589 : : do jatom=1,natom
1590 : : do katom=1,natom
1591 : : if (abs(asr4(jatom,katom,ii,jj,kk,ll)).gt.tol8) then
1592 : : write(std_out,'(a,1x,7(i3,1x),1(e17.10,1x))') '>>>>> WARNING (ASR4) --->',&
1593 : : & ii,jj,kk,ll,iatom,jatom,katom,asr4(jatom,katom,ii,jj,kk,ll)
1594 : : if (abs(asr4(jatom,katom,ii,jj,kk,ll)).gt.tol6)&
1595 : : & ABI_WARNING('The acoustic sum rule is not fulfilled (order 4)')
1596 : : end if
1597 : : end do !katom
1598 : : end do !jatom
1599 : : end do !ll
1600 : : end do !kk
1601 : : end do !jj
1602 : : end do !ii
1603 : : !FB! Check the rotational invariance (third order)
1604 : : !FB do jatom=1,natom
1605 : : !FB do alpha=1,3
1606 : : !FB do beta=1,3
1607 : : !FB do gama=1,3
1608 : : !FB do lambda=1,3
1609 : : !FB if (abs(rot3(jatom,alpha,beta,gama,lambda)).gt.tol8) then
1610 : : !FB write(std_out,'(a,6(i3,1x),1(e17.10,1x))') &
1611 : : !FB &'>>>>> WARNING (ROT3) ---> iatom,jatom,alpha,beta,gama,lambda,norm =',&
1612 : : !FB& iatom,jatom,alpha,beta,gama,lambda,rot3(jatom,alpha,beta,gama,lambda)
1613 : : !FB if (abs(rot3(jatom,alpha,beta,gama,lambda)).gt.tol6)&
1614 : : !FB& ABI_WARNING('The invariance under arbitrary rotation is not fulfilled (order 3)')
1615 : : !FB end if
1616 : : !FB end do !lambda
1617 : : !FB end do !gama
1618 : : !FB end do !beta
1619 : : !FB end do !alpha
1620 : : !FB end do !jatom
1621 : : end do !iatom
1622 : : ABI_FREE(asr4)
1623 : : !FB ABI_FREE(rot3)
1624 : : end if !order=4
1625 : :
1626 : 46 : end subroutine tdep_check_constraints
1627 : : !!***
1628 : :
1629 : : !=====================================================================================================
1630 : :
1631 : : !!****f* m_tdep_model/tdep_calc_gruneisen
1632 : : !! NAME
1633 : : !! tdep_calc_gruneisen
1634 : : !!
1635 : : !! FUNCTION
1636 : : !! Compute the mode Gruneisen parameters (as 3x3 tensors) at a given q-point,
1637 : : !! using the third-order IFC and the phonon eigenvectors and frequencies:
1638 : : !! Gamma_{nu}^{gamma delta}(q) = - 1/(6 omega_{nu}(q)^2)
1639 : : !! sum_{i in unitcell} sum_{j,k} sum_{alpha,beta}
1640 : : !! e*_{i alpha}(nu,q) e_{j beta}(nu,q) Phi3_{ijk}^{alpha beta gamma}
1641 : : !! F^{gamma delta} r_{ik}^{delta} exp(i q.R_j) / sqrt(M_i M_j)
1642 : : !! where r_{ik} is the vector between atoms i and k, and F is the strain
1643 : : !! tensor (currently the identity, so that only the diagonal elements of the
1644 : : !! Gruneisen tensors are non-zero).
1645 : : !!
1646 : : !! INPUTS
1647 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
1648 : : !! Eigen2nd = Eigen object containing the phonon frequencies, eigenvectors
1649 : : !! and dynamical matrices.
1650 : : !! iqpt = Index of the q-point in Eigen2nd.
1651 : : !! Invar = Input object containing the input variables.
1652 : : !! Phi3_ref(3,3,3,nshell) = Third-order IFC for the reference atoms of each shell.
1653 : : !! qpt_cart(3) = Cartesian coordinates of the q-point.
1654 : : !! Rlatt_cart(3,natom_unitcell,natom) = Cartesian lattice vector between
1655 : : !! a unitcell reference atom and the
1656 : : !! ideal site of each atom.
1657 : : !! Shell3at = Shell object for third-order interactions.
1658 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1659 : : !!
1660 : : !! OUTPUT
1661 : : !! Gruneisen(3*natom_unitcell,3,3) = Mode Gruneisen tensors at this q-point.
1662 : : !!
1663 : : !! SOURCE
1664 : :
1665 : 4635 : subroutine tdep_calc_gruneisen(SC,Eigen2nd,Gruneisen,iqpt,Invar,Phi3_ref,qpt_cart,Rlatt_cart,Shell3at,Sym)
1666 : :
1667 : : type(Symmetries_type),intent(in) :: Sym
1668 : : type(atdep_dataset_type),intent(in) :: Invar
1669 : : type(tdep_Supercell_type),intent(in) :: SC
1670 : : type(Shell_type),intent(in) :: Shell3at
1671 : : type(Eigen_type),intent(in) :: Eigen2nd
1672 : : integer,intent(in) :: iqpt
1673 : : double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
1674 : : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
1675 : : double precision,intent(in) :: qpt_cart(3)
1676 : : double complex ,intent(out) :: Gruneisen(3*Invar%natom_unitcell,3,3)
1677 : :
1678 : : integer :: ii,jj,kk,iatcell,jatom,katom,itypat,jtypat,natom_unitcell,ll
1679 : : integer :: imode,nmode,jat_mod,jatcell,iatshell,ishell,isym,itrans
1680 : : double precision :: phase
1681 : : double precision :: F_strain(3,3)
1682 : 4635 : double precision, allocatable :: Phi3_333(:,:,:)
1683 : 4635 : double complex, allocatable :: eigen_prod(:,:,:,:,:),eigenvec_loc(:,:),omega(:)
1684 : :
1685 : : ! Define quantities
1686 : 4635 : natom_unitcell=Invar%natom_unitcell
1687 : 4635 : nmode=3*Invar%natom_unitcell
1688 [ + + ]: 18540 : do ii=1,3
1689 [ + + ]: 60255 : do jj=1,3
1690 [ + + ]: 55620 : if (ii.eq.jj) then
1691 : : !FB F_strain(ii,jj)= 1./dsqrt(6.d0)
1692 : : !FB F_strain(ii,jj)=1./dsqrt(15.d0)
1693 : 13905 : F_strain(ii,jj)= 1.d0
1694 : : !FB F_strain(ii,jj)= 0.8
1695 : : else
1696 : : !FB F_strain(ii,jj)= 1./dsqrt(6.d0)/2.d0
1697 : : !FB F_strain(ii,jj)=2./dsqrt(15.d0)
1698 : 27810 : F_strain(ii,jj)=0.d0
1699 : : !FB F_strain(ii,jj)= 0.3
1700 : : end if
1701 : : end do
1702 : : end do
1703 : :
1704 : : ! Calculation of the eigenvectors product (and phonon frequencies)
1705 [ - + + - : 1710558 : ABI_MALLOC(eigen_prod, (3,3,natom_unitcell,natom_unitcell,nmode)) ; eigen_prod(:,:,:,:,:)=czero
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + +
+ ]
1706 [ + - - + : 214461 : ABI_MALLOC(eigenvec_loc,(3*natom_unitcell,3*natom_unitcell)) ; eigenvec_loc(:,:) =czero
+ - - + +
- + - + -
- + - + +
+ + + ]
1707 [ + + ]: 13416 : do iatcell=1,natom_unitcell
1708 [ + + ]: 31743 : do jatcell=1,natom_unitcell
1709 [ + + ]: 82089 : do ii=1,3
1710 [ + + ]: 238251 : do jj=1,3
1711 : : eigenvec_loc((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
1712 : 219924 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
1713 : : end do !ii
1714 : : end do !jj
1715 : : end do !jatcell
1716 : : end do !iatcell
1717 [ - + + - : 40248 : ABI_MALLOC(omega,(3*natom_unitcell)); omega(:)=czero
+ - + - -
+ - + +
+ ]
1718 [ + + ]: 13416 : do iatcell=1,natom_unitcell
1719 [ + + ]: 31743 : do jatcell=1,natom_unitcell
1720 [ + + ]: 82089 : do ii=1,3
1721 [ + + ]: 238251 : do jj=1,3
1722 [ + + ]: 1331811 : do imode=1,nmode
1723 : : eigen_prod(ii,jj,iatcell,jatcell,imode)=conjg(eigenvec_loc((iatcell-1)*3+ii,imode))*&
1724 : 1111887 : & eigenvec_loc((jatcell-1)*3+jj,imode)
1725 : : omega(imode)=omega(imode)+eigen_prod(ii,jj,iatcell,jatcell,imode)*&
1726 : : & dcmplx(Eigen2nd%dynmat(1,ii,iatcell,jj,jatcell,iqpt),&
1727 : 1276830 : & Eigen2nd%dynmat(2,ii,iatcell,jj,jatcell,iqpt))
1728 : : !FB write(std_out,'(a,5(1x,i4),2(1x,f15.3))') 'imode,iatcell,jatcell,ii,jj=',imode,iatcell,jatcell,ii,jj
1729 : : !FB eigen_prod(ii,jj,iatcell,jatcell,imode)
1730 : : end do !imode
1731 : : end do !ii
1732 : : end do !jj
1733 : : end do !jatcell
1734 : : end do !iatcell
1735 : :
1736 [ + + ]: 30978 : do imode=1,nmode
1737 [ - + ]: 30978 : if (abs(aimag(omega(imode))).gt.tol8) then
1738 : 0 : write(std_out,'(a,1x,e15.8,1x,a,i4)') '>>>> WARNING : Imaginary part of the phonon frequency is not zero (',&
1739 : 0 : & aimag(omega(imode)),') for mode :',imode
1740 : : !FB stop
1741 : : end if
1742 : : end do
1743 : : !FB write(std_out,'(a,1(1x,i4),6(1x,f15.3))') 'iqpt,Grun=',iqpt,(dsqrt(real(omega(imode))),imode=1,6)
1744 [ - + ]: 4635 : ABI_FREE(omega)
1745 : :
1746 : :
1747 : : ! Calculation of the Gruneisen
1748 [ - + + + : 185400 : ABI_MALLOC(Phi3_333,(3,3,3)) ; Phi3_333(:,:,:)=0.d0
+ + + + ]
1749 [ + + + + : 297342 : Gruneisen(:,:,:)=czero
+ + ]
1750 [ + + ]: 19691 : do ishell=1,Shell3at%nshell
1751 [ + + ]: 52449 : do iatcell=1,Invar%natom_unitcell
1752 : 32758 : itypat=Invar%typat_unitcell(iatcell)
1753 [ + + ]: 389103 : do iatshell=1,Shell3at%neighbours(iatcell,ishell)%n_interactions
1754 : 341289 : jatom=Shell3at%neighbours(iatcell,ishell)%atomj_in_shell(iatshell)
1755 : 341289 : katom=Shell3at%neighbours(iatcell,ishell)%atomk_in_shell(iatshell)
1756 : 341289 : isym =Shell3at%neighbours(iatcell,ishell)%sym_in_shell(iatshell)
1757 : 341289 : itrans=Shell3at%neighbours(iatcell,ishell)%transpose_in_shell(iatshell)
1758 : 341289 : jat_mod=mod(jatom+Invar%natom_unitcell-1,Invar%natom_unitcell)+1
1759 : 341289 : jtypat=Invar%typat_unitcell(jat_mod)
1760 : 341289 : phase=0.d0
1761 [ + + ]: 1365156 : do jj=1,3
1762 : 1365156 : phase=phase+2*pi*Rlatt_cart(jj,iatcell,jatom)*qpt_cart(jj)
1763 : : end do
1764 : 341289 : call tdep_build_phi3_333(isym,Phi3_ref(:,:,:,ishell),Phi3_333,Sym,itrans)
1765 [ + + ]: 1397914 : do ii=1,3
1766 [ + + ]: 4436757 : do jj=1,3
1767 [ + + ]: 13310271 : do kk=1,3
1768 [ + + ]: 39930813 : do ll=1,3
1769 [ + + ]: 226679877 : do imode=1,nmode
1770 : : Gruneisen(imode,kk,ll)=Gruneisen(imode,kk,ll)-dcmplx(Phi3_333(ii,jj,kk),0.d0)*&
1771 : : & eigen_prod(ii,jj,iatcell,jat_mod,imode)*&
1772 : : & exp(dcmplx(0.d0,phase))*dcmplx(F_strain(kk,ll),0.d0)*dcmplx(SC%dij(iatcell,katom,ll+1),0.d0)/&
1773 : : & dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)/&
1774 : 217465074 : & dcmplx(6*Eigen2nd%eigenval(imode,iqpt)**2,0.d0)
1775 : : !FB do ll=1,3
1776 : : !FB Gruneisen(imode,kk,ll)=Gruneisen(imode,kk,ll)-dcmplx(Phi3_333(ii,jj,kk),0.d0)*&
1777 : : !FB& eigen_prod(ii,jj,iatcell,jat_mod,imode)*&
1778 : : !FB& exp(dcmplx(0.d0,phase))*F_strain(kk,ll)*dcmplx(SC%dij(iatcell,katom,ll+1),0.d0)/&
1779 : : !FB& dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)/&
1780 : : !FB& dcmplx(6*Eigen2nd%eigenval(imode,iqpt)**2,0.d0)
1781 : : !FB!FB write(std_out,'(a,4(1x,i4),2(1x,f15.3))') 'ii,jj,kk,imode,Grun=',ii,jj,kk,imode,Gruneisen(imode,kk,kk)
1782 : : !FB end do !ll
1783 : : end do !imode
1784 : : end do !ll
1785 : : end do !kk
1786 : : end do !jj
1787 : : end do !ii
1788 : : end do !iatshell
1789 : : end do !iatcell
1790 : : end do !ishell
1791 : :
1792 : : !FB Gruneisen(:,:,:)=Gruneisen(:,:,:)*dcmplx(dsqrt(6.d0),0.d0)
1793 : : !FB do imode=1,nmode
1794 : : !FB write(std_out,'(a,1(1x,i4),6(1x,f15.3))') 'imode,Grun=',imode,(Gruneisen(imode,kk,kk),kk=1,3)
1795 : : !FB end do
1796 [ - + ]: 4635 : ABI_FREE(Phi3_333)
1797 [ - + ]: 4635 : ABI_FREE(eigen_prod)
1798 [ - + ]: 4635 : ABI_FREE(eigenvec_loc)
1799 : :
1800 : 4635 : end subroutine tdep_calc_gruneisen
1801 : : !!***
1802 : :
1803 : : !=====================================================================================================
1804 : :
1805 : : !!****f* m_tdep_model/tdep_calc_alpha_gamma
1806 : : !! NAME
1807 : : !! tdep_calc_alpha_gamma
1808 : : !!
1809 : : !! FUNCTION
1810 : : !! Compute the thermodynamic Gruneisen parameter (Gamma), the thermal
1811 : : !! expansion (alpha_v), the thermal energy and the thermal pressure (P_th),
1812 : : !! as well as the heat capacities at constant volume and pressure and the
1813 : : !! isentropic bulk modulus, by integrating the mode Gruneisen parameters
1814 : : !! (see tdep_calc_gruneisen) over the irreducible Brillouin zone.
1815 : : !! These quantities are computed at Invar%temperature and written in the
1816 : : !! output file. They are also computed, under the harmonic approximation,
1817 : : !! as a function of temperature (from 10 K to 10000 K) and written in the
1818 : : !! _thermo3.dat and _alpha_gamma.dat files.
1819 : : !!
1820 : : !! INPUTS
1821 : : !! Invar = Input object containing the input variables.
1822 : : !! MPIdata = Info on MPI parallelism.
1823 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
1824 : : !! Eigen2nd = Eigen object containing the phonon frequencies and eigenvectors
1825 : : !! on the irreducible q-points of Qbz.
1826 : : !! Qbz = Qbz object containing the q-points of the irreducible Brillouin zone
1827 : : !! and their weights.
1828 : : !! Rlatt_cart(3,natom_unitcell,natom) = Cartesian lattice vector between
1829 : : !! a unitcell reference atom and the
1830 : : !! ideal site of each atom.
1831 : : !! Shell3at = Shell object for third-order interactions.
1832 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1833 : : !!
1834 : : !! OUTPUT
1835 : : !!
1836 : : !! SIDE EFFECTS
1837 : : !! Model = TDEP Model object. On entry, Model%Phi3, Model%Sij, Model%ucvol,
1838 : : !! Model%BulkModulus_T, Model%Shear and Model%Density must already
1839 : : !! be set (see tdep_model_calc_elastic). On exit,
1840 : : !! Model%BulkModulus_S and Model%HeatCapa_P are computed and stored.
1841 : : !!
1842 : : !! SOURCE
1843 : :
1844 : 16 : subroutine tdep_calc_alpha_gamma(Model,Invar,MPIdata,SC,Eigen2nd,Qbz,Rlatt_cart,Shell3at,Sym)
1845 : :
1846 : : class(tdep_Model_type),intent(inout) :: Model
1847 : : type(atdep_dataset_type),intent(in) :: Invar
1848 : : type(MPI_enreg_type), intent(in) :: MPIdata
1849 : : type(tdep_Supercell_type),intent(in) :: SC
1850 : : type(Eigen_type),intent(in) :: Eigen2nd
1851 : : type(Shell_type),intent(in) :: Shell3at
1852 : : type(Symmetries_type),intent(in) :: Sym
1853 : : type(Qbz_type),intent(in) :: Qbz
1854 : : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
1855 : :
1856 : : integer :: iq_ibz,ii,jj,nmode,ntemp,itemp,kk,imode,ll
1857 : : integer :: alpha,beta
1858 : : double precision :: k_B,wovert,xx,C_v,E_th,Vp,Vs,P_th2
1859 : 16 : double precision, allocatable :: heatcapa(:),grun_thermo(:,:,:),u_vib(:)
1860 : 16 : double precision, allocatable :: heatcapa_HA(:,:),grun_thermo_HA(:,:,:,:),p_thermo1(:,:,:)
1861 : 16 : double precision, allocatable :: p_thermo_HA(:,:,:,:),u_vib_HA(:,:)
1862 : 16 : double precision, allocatable :: p_thermo2(:),tmp(:)
1863 : : double precision :: Gama_tensor(3,3),alpha_v_tensor(3,3),P_th_tensor(3,3),P_th1_tensor(3,3)
1864 : 16 : double complex, allocatable :: Gruneisen(:,:,:),Grun_mean(:)
1865 : : character(len=500) :: message
1866 : : integer :: ierr
1867 : :
1868 : 16 : ierr = 0
1869 : :
1870 [ + - - + : 1152 : ABI_MALLOC(Gruneisen,(3*Invar%natom_unitcell,3,3)); Gruneisen(:,:,:)=czero
+ - + - +
- + - + -
- + - + +
+ + + +
+ ]
1871 [ - + + - : 144 : ABI_MALLOC(Grun_mean,(3*Invar%natom_unitcell)) ; Grun_mean(:) =czero
+ - + - -
+ - + +
+ ]
1872 : :
1873 : 16 : ntemp=1000
1874 : 16 : nmode=3*Invar%natom_unitcell
1875 : 16 : k_B=kb_HaK*Ha_eV
1876 [ - + + - : 144 : ABI_MALLOC(heatcapa,(nmode)) ; heatcapa (:) =0.d0
+ - + - -
+ - + +
+ ]
1877 [ - + + - : 1120 : ABI_MALLOC(grun_thermo,(nmode,3,3)); grun_thermo(:,:,:)=0.d0
+ - + - +
- - + - +
+ + + + +
+ ]
1878 [ - + + - : 1088 : ABI_MALLOC(p_thermo1,(nmode,3,3)) ; p_thermo1 (:,:,:)=0.d0
- + - + +
+ + + +
+ ]
1879 [ - + + + ]: 16016 : ABI_MALLOC(p_thermo2,(ntemp)) ; p_thermo2 (:) =0.d0
1880 [ - + + - : 128 : ABI_MALLOC(u_vib,(nmode)) ; u_vib (:) =0.d0
- + - + +
+ ]
1881 [ - + + - : 112064 : ABI_MALLOC(heatcapa_HA,(nmode,ntemp)) ; heatcapa_HA (:,:) =0.d0
+ - + - +
- - + - +
+ + + + ]
1882 [ - + + - : 1008288 : ABI_MALLOC(grun_thermo_HA,(nmode,ntemp,3,3)); grun_thermo_HA(:,:,:,:)=0.d0
+ - + - +
- + - + -
- + - + +
+ + + + +
+ + ]
1883 [ - + + - : 1008240 : ABI_MALLOC(p_thermo_HA,(nmode,ntemp,3,3)) ; p_thermo_HA (:,:,:,:)=0.d0
+ - - + -
+ + + + +
+ + + + ]
1884 [ - + + - : 112048 : ABI_MALLOC(u_vib_HA,(nmode,ntemp)) ; u_vib_HA (:,:) =0.d0
+ - - + -
+ + + +
+ ]
1885 : : ! Loop over irreducible q-points
1886 : : ! =======================
1887 [ + + ]: 1568 : do iq_ibz=1,Qbz%nqibz
1888 : :
1889 : : ! Compute the gruneisen for this q-point
1890 : : ! ======================================
1891 [ + + + + ]: 6208 : if ((sum(abs(Qbz%qibz_cart(:,iq_ibz)))).lt.tol8) cycle ! G point
1892 : : call tdep_calc_gruneisen(SC,Eigen2nd,Gruneisen,iq_ibz,Invar,&
1893 : 1536 : & Model%Phi3,Qbz%qibz_cart(:,iq_ibz),Rlatt_cart,Shell3at,Sym)
1894 [ + + ]: 10572 : Grun_mean(:)=czero
1895 [ + + ]: 10572 : do ii=1,3*Invar%natom_unitcell
1896 [ + + ]: 37680 : do jj=1,3
1897 [ + + ]: 117468 : do kk=1,3
1898 : 108432 : Grun_mean(ii)=Grun_mean(ii)+Gruneisen(ii,jj,kk)
1899 : : end do
1900 : : end do
1901 : : end do
1902 [ + + - + ]: 10572 : if (sum(abs(aimag(Grun_mean(:)))).gt.tol8) then
1903 [ # # # # ]: 0 : write(message,'(i5,1x,100(e15.6,1x))') iq_ibz,( real(Grun_mean(ii)),ii=1,nmode)
1904 : 0 : ABI_ERROR_NOSTOP(message,ierr)
1905 [ # # # # ]: 0 : write(message,'(i5,1x,100(e15.6,1x))') iq_ibz,(aimag(Grun_mean(ii)),ii=1,nmode)
1906 : 0 : ABI_ERROR_NOSTOP(message,ierr)
1907 : 0 : ABI_ERROR('tdep_calc_alpha_gamma : The imaginary part of the Gruneisen is not equal to zero')
1908 : : end if
1909 : : ! If the Gruneisen is isotropic
1910 : : !FB do ii=1,3
1911 : : !FB tmp1=0.d0
1912 : : !FB do jj=1,3
1913 : : !FB tmp1=tmp1 + real(Gruneisen(ii,jj,jj))
1914 : : !FB end do
1915 : : !FB do jj=1,3
1916 : : !FB Gruneisen(ii,jj,jj)=dcmplx(tmp1/3.d0,0.d0)
1917 : : !FB end do
1918 : : !FB write(std_out,'(2(i5,1x),100(e15.6,1x))') iq_ibz,ii,Gruneisen(ii,jj,jj)
1919 : : !FB end do
1920 : :
1921 : : ! Compute the heat capacity and thermodynamical gruneisen parameter at present temperature
1922 : : ! ========================================================================================
1923 : 1536 : wovert=1.d0/(2*Invar%temperature*k_B)
1924 [ + + ]: 10572 : do ii=1,nmode
1925 : 9036 : xx=Eigen2nd%eigenval(ii,iq_ibz)*Ha_eV
1926 [ - + ]: 9036 : if (xx.le.0) cycle
1927 : 9036 : heatcapa(ii) =heatcapa(ii) +(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)
1928 : : ! Gamma= sum_i Gamma_i*C_Vi / C_V
1929 [ + + + + ]: 117468 : grun_thermo(ii,:,:)=grun_thermo(ii,:,:)+(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)*real(Gruneisen(ii,:,:))*3.d0
1930 : : ! P = sum_i Gamma_i*E_i / V
1931 [ + + + + ]: 117468 : p_thermo1(ii,:,:) =p_thermo1(ii,:,:) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)*real(Gruneisen(ii,:,:))*3.d0
1932 : 10572 : u_vib(ii) =u_vib(ii) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)
1933 : : end do
1934 : : ! Compute the heat capacity and thermodynamical gruneisen parameter as a function of temperature
1935 : : ! ==============================================================================================
1936 [ + + ]: 1537552 : do itemp=1,ntemp
1937 : 1536000 : wovert=1.d0/(2*real(itemp)*10*k_B)
1938 [ + + ]: 10573552 : do ii=1,nmode
1939 : 9036000 : xx=Eigen2nd%eigenval(ii,iq_ibz)*Ha_eV
1940 [ - + ]: 9036000 : if (xx.le.0) cycle
1941 : 9036000 : heatcapa_HA(ii,itemp) =heatcapa_HA(ii,itemp) +(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)
1942 : : grun_thermo_HA(ii,itemp,:,:)=grun_thermo_HA(ii,itemp,:,:)+(wovert*xx/sinh(wovert*xx))**2*Qbz%wtqibz(iq_ibz)&
1943 [ + + + + ]: 117468000 : & *real(Gruneisen(ii,:,:))*3.d0
1944 : : p_thermo_HA(ii,itemp,:,:) =p_thermo_HA(ii,itemp,:,:) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)&
1945 [ + + + + ]: 117468000 : & *real(Gruneisen(ii,:,:))*3.d0
1946 : 10572000 : u_vib_HA(ii,itemp) =u_vib_HA(ii,itemp) + (xx/2.d0/tanh(wovert*xx)) *Qbz%wtqibz(iq_ibz)
1947 : : end do
1948 : : end do
1949 : : end do
1950 [ - + ]: 16 : ABI_FREE(Gruneisen)
1951 [ - + ]: 16 : ABI_FREE(Grun_mean)
1952 : : ! Compute the pressure as the integral of Gamma*C_v overt T
1953 [ - + + + ]: 16016 : ABI_MALLOC(tmp,(ntemp)) ; tmp(:)=0.d0
1954 [ + + ]: 16016 : do itemp=1,ntemp
1955 [ + + ]: 64016 : do ii=1,3
1956 [ + + ]: 208000 : do jj=1,3
1957 [ + + ]: 1056000 : tmp(itemp)=tmp(itemp)+sum(grun_thermo_HA(:,itemp,ii,jj))/3.d0
1958 : : end do
1959 : : end do
1960 : : end do
1961 : 16 : call simpson_int(ntemp,10.d0,tmp,p_thermo2)
1962 [ - + ]: 16 : ABI_FREE(tmp)
1963 : :
1964 [ + - ]: 16 : if (MPIdata%iam_master) then
1965 : 16 : open(unit=20,file=trim(Invar%output_prefix)//'_thermo3.dat')
1966 : 16 : open(unit=21,file=trim(Invar%output_prefix)//'_alpha_gamma.dat')
1967 : 16 : write(20,'(a)')'# T(K) C_v(k_B/fu) Gamma alpha_v*10^6(K^-1) E_th(eV) P_th(GPa)'
1968 : 16 : write(20,'(a,72x,a)')'#',' ----------------------------------------------'
1969 : 16 : write(20,'(a,72x,a)')'#',' {sum G_i.U_iV} {int G.C_v/V dT} {G.U/V}'
1970 : 16 : write(21,'(a)')'# T(K) Gamma_11 Gamma_22 Gamma_33 alpha_11 alpha_22 alpha_33',&
1971 : 32 : & ' alpha_12 alpha_13 alpha_23'
1972 : 16 : write(20,'(2a)')'# --------------------------------------------------------------------------------------------------',&
1973 : 32 : '--------------------'
1974 : : end if
1975 [ + + ]: 16016 : do itemp=1,ntemp
1976 [ + + ]: 112000 : C_v =sum(heatcapa_HA (:,itemp))
1977 [ + + ]: 112000 : E_th =sum(u_vib_HA (:,itemp))
1978 : 16000 : alpha_v_tensor(:,:)=zero
1979 : 16000 : Gama_tensor (:,:)=zero
1980 : 16000 : P_th1_tensor (:,:)=zero
1981 : 16000 : P_th_tensor (:,:)=zero
1982 [ + + ]: 112000 : do imode=1,nmode
1983 [ + + + + ]: 1248000 : Gama_tensor (:,:)=Gama_tensor (:,:)+grun_thermo_HA(imode,itemp,:,:)/C_v
1984 [ + + + + ]: 1264000 : P_th1_tensor(:,:)=P_th1_tensor(:,:)+p_thermo_HA (imode,itemp,:,:)/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
1985 : : end do
1986 : : ! alpha=sum{I=1,6;J=1,6} S(I,J) Gama_tensor (J)
1987 [ + + ]: 64000 : do ii=1,3
1988 [ + + ]: 208000 : do jj=1,3
1989 [ + + ]: 144000 : if (ii.eq.1.and.jj.eq.1) alpha=1
1990 [ + + ]: 144000 : if (ii.eq.2.and.jj.eq.2) alpha=2
1991 [ + + ]: 144000 : if (ii.eq.3.and.jj.eq.3) alpha=3
1992 [ + + ]: 144000 : if (ii.eq.2.and.jj.eq.3) alpha=4
1993 [ + + ]: 144000 : if (ii.eq.1.and.jj.eq.3) alpha=5
1994 [ + + ]: 144000 : if (ii.eq.1.and.jj.eq.2) alpha=6
1995 [ + + ]: 144000 : if (ii.eq.3.and.jj.eq.2) alpha=4
1996 [ + + ]: 144000 : if (ii.eq.3.and.jj.eq.1) alpha=5
1997 [ + + ]: 144000 : if (ii.eq.2.and.jj.eq.1) alpha=6
1998 [ + + ]: 576000 : do kk=1,3
1999 [ + + ]: 1872000 : do ll=1,3
2000 [ + + ]: 1296000 : if (kk.eq.1.and.ll.eq.1) beta=1
2001 [ + + ]: 1296000 : if (kk.eq.2.and.ll.eq.2) beta=2
2002 [ + + ]: 1296000 : if (kk.eq.3.and.ll.eq.3) beta=3
2003 [ + + ]: 1296000 : if (kk.eq.2.and.ll.eq.3) beta=4
2004 [ + + ]: 1296000 : if (kk.eq.1.and.ll.eq.3) beta=5
2005 [ + + ]: 1296000 : if (kk.eq.1.and.ll.eq.2) beta=6
2006 [ + + ]: 1296000 : if (kk.eq.3.and.ll.eq.2) beta=4
2007 [ + + ]: 1296000 : if (kk.eq.3.and.ll.eq.1) beta=5
2008 [ + + ]: 1296000 : if (kk.eq.2.and.ll.eq.1) beta=6
2009 : : alpha_v_tensor(ii,jj)=alpha_v_tensor(ii,jj)+Model%Sij(alpha,beta)*Gama_tensor(kk,ll)*C_v*kb_HaK*Ha_J/1.d9/&
2010 : 1728000 : & (Model%ucvol*Bohr_Ang**3*1.d-30)
2011 : : end do
2012 : : end do
2013 : 192000 : P_th_tensor(ii,jj) =Gama_tensor(ii,jj)*E_th/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
2014 : : end do
2015 : : end do
2016 : 16000 : P_th2 = p_thermo2(itemp)*k_B/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
2017 [ + - ]: 16016 : if (MPIdata%iam_master) then
2018 : 16000 : 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,&
2019 : 16000 : & (alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*1.d6,&
2020 : 16000 : & E_th,(P_th1_tensor(1,1)+P_th1_tensor(2,2)+P_th1_tensor(3,3))/3.d0,&
2021 : 32000 : & P_th2,(P_th_tensor(1,1)+P_th_tensor(2,2)+P_th_tensor(3,3))/3.d0
2022 : 16000 : write(21,'(1x,i5,9(1x,f15.3))') itemp*10,Gama_tensor(1,1),Gama_tensor(2,2),Gama_tensor(3,3),&
2023 : 16000 : & alpha_v_tensor(1,1)*1.d6,alpha_v_tensor(2,2)*1.d6,alpha_v_tensor(3,3)*1.d6,&
2024 : 32000 : & alpha_v_tensor(1,2)*1.d6,alpha_v_tensor(1,3)*1.d6,alpha_v_tensor(2,3)*1.d6
2025 : : end if
2026 : : end do
2027 [ + - ]: 16 : if (MPIdata%iam_master) then
2028 : 16 : close(20)
2029 : 16 : close(21)
2030 : : end if
2031 [ - + ]: 16 : ABI_FREE(heatcapa_HA)
2032 [ - + ]: 16 : ABI_FREE(grun_thermo_HA)
2033 [ - + ]: 16 : ABI_FREE(p_thermo_HA)
2034 [ - + ]: 16 : ABI_FREE(u_vib_HA)
2035 : :
2036 : : ! Summary
2037 : : ! =======
2038 : 16 : write(ab_out,*) ' '
2039 : 16 : write(ab_out,*) '#############################################################################'
2040 : 16 : write(ab_out,*) '####### Gruneisen parameter, Thermal expansion, Thermal pressure... #########'
2041 : 16 : write(ab_out,*) '#############################################################################'
2042 [ + + ]: 112 : C_v =sum(heatcapa(:))
2043 [ + + ]: 112 : E_th =sum(u_vib (:))
2044 : 16 : alpha_v_tensor(:,:)=zero
2045 : 16 : Gama_tensor (:,:)=zero
2046 : 16 : P_th1_tensor (:,:)=zero
2047 : 16 : P_th_tensor (:,:)=zero
2048 [ + + ]: 112 : do imode=1,nmode
2049 [ + + + + ]: 1248 : Gama_tensor (:,:)=Gama_tensor (:,:)+grun_thermo(imode,:,:)/C_v
2050 [ + + + + ]: 1264 : P_th1_tensor(:,:)=P_th1_tensor(:,:)+p_thermo1 (imode,:,:)/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
2051 : : end do
2052 : : ! alpha=sum{I=1,6;J=1,6} S(I,J) Gama_tensor (J)
2053 [ + + ]: 64 : do ii=1,3
2054 [ + + ]: 208 : do jj=1,3
2055 [ + + ]: 144 : if (ii.eq.1.and.jj.eq.1) alpha=1
2056 [ + + ]: 144 : if (ii.eq.2.and.jj.eq.2) alpha=2
2057 [ + + ]: 144 : if (ii.eq.3.and.jj.eq.3) alpha=3
2058 [ + + ]: 144 : if (ii.eq.2.and.jj.eq.3) alpha=4
2059 [ + + ]: 144 : if (ii.eq.1.and.jj.eq.3) alpha=5
2060 [ + + ]: 144 : if (ii.eq.1.and.jj.eq.2) alpha=6
2061 [ + + ]: 144 : if (ii.eq.3.and.jj.eq.2) alpha=4
2062 [ + + ]: 144 : if (ii.eq.3.and.jj.eq.1) alpha=5
2063 [ + + ]: 144 : if (ii.eq.2.and.jj.eq.1) alpha=6
2064 [ + + ]: 576 : do kk=1,3
2065 [ + + ]: 1872 : do ll=1,3
2066 [ + + ]: 1296 : if (kk.eq.1.and.ll.eq.1) beta=1
2067 [ + + ]: 1296 : if (kk.eq.2.and.ll.eq.2) beta=2
2068 [ + + ]: 1296 : if (kk.eq.3.and.ll.eq.3) beta=3
2069 [ + + ]: 1296 : if (kk.eq.2.and.ll.eq.3) beta=4
2070 [ + + ]: 1296 : if (kk.eq.1.and.ll.eq.3) beta=5
2071 [ + + ]: 1296 : if (kk.eq.1.and.ll.eq.2) beta=6
2072 [ + + ]: 1296 : if (kk.eq.3.and.ll.eq.2) beta=4
2073 [ + + ]: 1296 : if (kk.eq.3.and.ll.eq.1) beta=5
2074 [ + + ]: 1296 : if (kk.eq.2.and.ll.eq.1) beta=6
2075 : : alpha_v_tensor(ii,jj)=alpha_v_tensor(ii,jj)+Model%Sij(alpha,beta)*Gama_tensor(kk,ll)*C_v*kb_HaK*Ha_J/1.d9/&
2076 : 1728 : & (Model%ucvol*Bohr_Ang**3*1.d-30)
2077 : : end do
2078 : : end do
2079 : 192 : P_th_tensor(ii,jj)=Gama_tensor(ii,jj)*E_th/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
2080 : : end do
2081 : : end do
2082 : 16 : P_th2 =p_thermo2(int(Invar%temperature/10))*k_B/(Model%ucvol*Bohr_Ang**3*1.d-30)*e_Cb/10**9
2083 : : Model%BulkModulus_S=Model%BulkModulus_T*(1.+(alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*&
2084 : 16 : & (Gama_tensor(1,1)+ Gama_tensor(2,2)+ Gama_tensor(3,3))/3.d0*Invar%temperature)
2085 : 16 : Model%HeatCapa_P=C_v*Model%BulkModulus_S/Model%BulkModulus_T
2086 : 16 : Vp=dsqrt(1.d9*(Model%BulkModulus_S+4.d0*Model%Shear/3.d0)/Model%Density)
2087 : 16 : Vs=dsqrt(1.d9*Model%Shear/Model%Density)
2088 : 16 : write(ab_out,'(a)') ' See the gruneisen.dat, alpha_gamma.dat and thermo3.dat files'
2089 : 16 : write(ab_out,'(a)') ' '
2090 : 16 : write(ab_out,'(a,1x,f15.3)') ' Gruneisen parameter : Gamma=',&
2091 : 32 : & (Gama_tensor(1,1)+Gama_tensor(2,2)+Gama_tensor(3,3))/3.d0
2092 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' The Gruneisen matrix is |',&
2093 : 32 : & Gama_tensor(1,1),Gama_tensor(1,2),Gama_tensor(1,3),'|'
2094 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2095 : 32 : & Gama_tensor(2,1),Gama_tensor(2,2),Gama_tensor(2,3),'|'
2096 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2097 : 32 : & Gama_tensor(3,1),Gama_tensor(3,2),Gama_tensor(3,3),'|'
2098 : 16 : write(ab_out,'(a,1x,f15.3)') ' Thermal expansion (K^{-1}*10^6) : alpha_v=',&
2099 : 32 : & (alpha_v_tensor(1,1)+alpha_v_tensor(2,2)+alpha_v_tensor(3,3))*1.d6
2100 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' The thermal expansion matrix is |',&
2101 : 32 : & alpha_v_tensor(1,1)*1.d6,alpha_v_tensor(1,2)*1.d6,alpha_v_tensor(1,3)*1.d6,'|'
2102 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2103 : 32 : & alpha_v_tensor(2,1)*1.d6,alpha_v_tensor(2,2)*1.d6,alpha_v_tensor(2,3)*1.d6,'|'
2104 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2105 : 32 : & alpha_v_tensor(3,1)*1.d6,alpha_v_tensor(3,2)*1.d6,alpha_v_tensor(3,3)*1.d6,'|'
2106 : 16 : write(ab_out,'(a)') ' Thermal pressure (in GPa) : '
2107 : 16 : write(ab_out,'(a)') ' ------- w intrinsic effects and w ZPE --------'
2108 : 16 : write(ab_out,'(a,1x,f15.3)') ' P_th=sum_i Gamma_i*E_i/V =',&
2109 : 32 : & (P_th1_tensor(1,1)+P_th1_tensor(2,2)+P_th1_tensor(3,3))/3.d0
2110 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' The thermal pressure matrix is|',&
2111 : 32 : & P_th1_tensor(1,1),P_th1_tensor(1,2),P_th1_tensor(1,3),'|'
2112 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2113 : 32 : & P_th1_tensor(2,1),P_th1_tensor(2,2),P_th1_tensor(2,3),'|'
2114 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2115 : 32 : & P_th1_tensor(3,1),P_th1_tensor(3,2),P_th1_tensor(3,3),'|'
2116 : 16 : write(ab_out,'(a)') ' ------- w intrinsic effects and w/o ZPE --------'
2117 : 16 : write(ab_out,'(a,1x,f15.3)') ' P_th=integ{Gamma*C_v/V dT} =',P_th2
2118 : 16 : write(ab_out,'(a)') ' ------- w/o intrinsic effects and w/o ZPE --------'
2119 : 16 : write(ab_out,'(a,1x,f15.3)') ' P_th=Gamma*E_th/V =',&
2120 : 32 : & (P_th_tensor(1,1)+P_th_tensor(2,2)+P_th_tensor(3,3))/3.d0
2121 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' The thermal pressure matrix is|',&
2122 : 32 : & P_th_tensor(1,1),P_th_tensor(1,2),P_th_tensor(1,3),'|'
2123 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2124 : 32 : & P_th_tensor(2,1),P_th_tensor(2,2),P_th_tensor(2,3),'|'
2125 : 16 : write(ab_out,'(a,3(1x,f8.3),a)') ' |',&
2126 : 32 : & P_th_tensor(3,1),P_th_tensor(3,2),P_th_tensor(3,3),'|'
2127 : 16 : write(ab_out,'(a,1x,f15.3)') ' Volume (bohr^3 per unit cell): V=',Model%ucvol
2128 : 16 : write(ab_out,'(a,1x,f15.3)') ' Thermal energy (eV) : E_th=',E_th
2129 : 16 : write(ab_out,'(a,1x,f15.3)') ' Heat capacity at constant V (k_B/f.u.): C_v=',C_v
2130 : 16 : write(ab_out,'(a,1x,f15.3)') ' Heat capacity at constant P (k_B/f.u.): C_p=',Model%HeatCapa_P
2131 : 16 : write(ab_out,'(a,1x,f15.3)') ' Isothermal Bulk Modulus (GPa): B_T=',Model%BulkModulus_T
2132 : 16 : write(ab_out,'(a,1x,f15.3)') ' Isentropic Bulk Modulus (GPa): B_S=',Model%BulkModulus_S
2133 : 16 : write(ab_out,'(a,1x,f15.3)') ' Longitudinal sound velocity (m.s-1): Vp=',Vp
2134 : 16 : write(ab_out,'(a,1x,f15.3)') ' Transverse sound velocity (m.s-1): Vs=',Vs
2135 [ - + ]: 16 : ABI_FREE(heatcapa)
2136 [ - + ]: 16 : ABI_FREE(grun_thermo)
2137 [ - + ]: 16 : ABI_FREE(p_thermo1)
2138 : 16 : ABI_FREE(p_thermo2)
2139 [ - + ]: 16 : ABI_FREE(u_vib)
2140 : :
2141 : 16 : end subroutine tdep_calc_alpha_gamma
2142 : : !!***
2143 : :
2144 : : !=====================================================================================================
2145 : :
2146 : : !!****f* m_tdep_model/tdep_calc_write_gruneisen
2147 : : !! NAME
2148 : : !! tdep_calc_write_gruneisen
2149 : : !!
2150 : : !! FUNCTION
2151 : : !! Compute the mode Gruneisen parameters (see tdep_calc_gruneisen) along the
2152 : : !! q-point path, and write them in the _gruneisen.dat file (sum of the
2153 : : !! components of the Gruneisen tensor for each mode) and in the
2154 : : !! _gruneisen-ij.dat file (all
2155 : : !! components of the Gruneisen tensor for each mode).
2156 : : !! The Gamma point (and its periodic images) is skipped.
2157 : : !! Only the master process performs this task.
2158 : : !!
2159 : : !! INPUTS
2160 : : !! Model = TDEP Model object, containing the third-order IFC (Model%Phi3).
2161 : : !! Invar = Input object containing the input variables.
2162 : : !! MPIdata = Info on MPI parallelism.
2163 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
2164 : : !! Eigen2nd = Eigen object containing the phonon frequencies and eigenvectors
2165 : : !! on the q-points of Qpt.
2166 : : !! Qpt = Qpoints object containing the q-point path.
2167 : : !! Rlatt_cart(3,natom_unitcell,natom) = Cartesian lattice vector between
2168 : : !! a unitcell reference atom and the
2169 : : !! ideal site of each atom.
2170 : : !! Shell3at = Shell object for third-order interactions.
2171 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
2172 : : !!
2173 : : !! OUTPUT
2174 : : !! Only writing.
2175 : : !!
2176 : : !! SOURCE
2177 : :
2178 : 16 : subroutine tdep_calc_write_gruneisen(Model,Invar,MPIdata,SC,Eigen2nd,Qpt,Rlatt_cart,Shell3at,Sym)
2179 : :
2180 : : class(tdep_Model_type),intent(in) :: Model
2181 : : type(atdep_dataset_type),intent(in) :: Invar
2182 : : type(MPI_enreg_type), intent(in) :: MPIdata
2183 : : type(tdep_Supercell_type),intent(in) :: SC
2184 : : type(Eigen_type),intent(in) :: Eigen2nd
2185 : : type(Qpoints_type),intent(in) :: Qpt
2186 : : type(Shell_type),intent(in) :: Shell3at
2187 : : type(Symmetries_type),intent(in) :: Sym
2188 : : real(dp),intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
2189 : :
2190 : : integer :: iqpt,nmode,ii,jj,kk
2191 : : double precision :: qpt_cart(3)
2192 : 16 : double complex, allocatable :: Gruneisen(:,:,:)
2193 : 16 : double complex, allocatable :: Grun_mean(:)
2194 : : character(len=500) :: message
2195 : : integer :: ierr
2196 : :
2197 [ - + ]: 16 : if (.not. MPIdata%iam_master) return
2198 : :
2199 : 16 : ierr = 0
2200 : :
2201 : 16 : nmode=3*Invar%natom_unitcell
2202 [ + - - + : 1152 : ABI_MALLOC(Gruneisen,(3*Invar%natom_unitcell,3,3)); Gruneisen(:,:,:)=czero
+ - + - +
- + - + -
- + - + +
+ + + +
+ ]
2203 [ - + + - : 144 : ABI_MALLOC(Grun_mean,(3*Invar%natom_unitcell)) ; Grun_mean(:) =czero
+ - + - -
+ - + +
+ ]
2204 : 16 : open(unit=53,file=trim(Invar%output_prefix)//'_gruneisen.dat')
2205 : 16 : open(unit=54,file=trim(Invar%output_prefix)//'_gruneisen-ij.dat')
2206 : :
2207 [ + + ]: 3147 : do iqpt=1,Qpt%nqpt
2208 : :
2209 : : ! Compute the Gruneisen
2210 : : ! =====================
2211 [ + + ]: 20630 : Grun_mean(:)=czero
2212 [ + + ]: 12524 : qpt_cart(:)=Qpt%qpt_cart(:,iqpt)
2213 [ + + + + ]: 12524 : if ((sum(abs(Qpt%qpt_red(:,iqpt)))).lt.tol8) cycle ! G point
2214 [ + + + + ]: 12404 : if (abs(sum(Qpt%qpt_red(:,iqpt)**2)-1.d0).lt.tol8) cycle ! Gp point
2215 : 3099 : call tdep_calc_gruneisen(SC,Eigen2nd,Gruneisen,iqpt,Invar,Model%Phi3,qpt_cart,Rlatt_cart,Shell3at,Sym)
2216 [ + + ]: 20406 : do ii=1,3*Invar%natom_unitcell
2217 [ + + ]: 72327 : do jj=1,3
2218 [ + + ]: 224991 : do kk=1,3
2219 : 155763 : write(54,'(4(i5,1x),500(e15.6,1x))') iqpt,ii,jj,kk,Gruneisen(ii,jj,kk)
2220 [ - + ]: 207684 : if (abs(aimag(Gruneisen(ii,jj,kk))).gt.tol8) then
2221 : 0 : ABI_WARNING('Real part :')
2222 : 0 : write(message,'(4(i5,1x),100(e15.6,1x))') iqpt,ii,jj,kk,real(Gruneisen(ii,jj,kk))
2223 : 0 : ABI_WARNING(message)
2224 : 0 : ABI_WARNING('Imaginary part :')
2225 : 0 : write(message,'(4(i5,1x),100(e15.6,1x))') iqpt,ii,jj,kk,aimag(Gruneisen(ii,jj,kk))
2226 : 0 : ABI_WARNING(message)
2227 : 0 : ABI_WARNING('tdep_calc_write_gruneisen : The imaginary part of the Gruneisen is not equal to zero')
2228 : : end if
2229 : : end do
2230 : : end do
2231 : : end do
2232 [ + + ]: 20406 : do ii=1,3*Invar%natom_unitcell
2233 [ + + ]: 72327 : do jj=1,3
2234 [ + + ]: 224991 : do kk=1,3
2235 : 207684 : Grun_mean(ii)=Grun_mean(ii)+Gruneisen(ii,jj,kk)
2236 : : end do
2237 : : end do
2238 : : end do
2239 : :
2240 : : ! Write the Gruneisen
2241 : : ! ===================
2242 [ + + - + ]: 20422 : if (sum(abs(aimag(Grun_mean(:)))).gt.3*Invar%natom_unitcell*tol8) then
2243 : 0 : ABI_WARNING('Real part :')
2244 [ # # # # ]: 0 : write(message,'(i5,1x,100(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode)
2245 : 0 : ABI_WARNING(message)
2246 : 0 : ABI_WARNING('Imaginary part :')
2247 [ # # # # ]: 0 : write(message,'(i5,1x,100(e15.6,1x))') iqpt,(aimag(Grun_mean(ii)),ii=1,nmode)
2248 : 0 : ABI_WARNING(message)
2249 : 0 : ABI_WARNING('tdep_calc_write_gruneisen : The imaginary part of Grun_mean is not equal to zero')
2250 : : else
2251 : : !FB write(53,'(i5,1x,500(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode),&
2252 : : !FB ((real(Grun_shell(ii,jj)),ii=1,nmode),jj=1,Shell3at%nshell)
2253 [ + - + + ]: 20406 : write(53,'(i5,1x,500(e15.6,1x))') iqpt,(real(Grun_mean(ii)),ii=1,nmode)
2254 : : end if
2255 : : end do
2256 : 16 : close(53)
2257 : 16 : close(54)
2258 [ - + ]: 16 : ABI_FREE(Grun_mean)
2259 [ - + ]: 16 : ABI_FREE(Gruneisen)
2260 : :
2261 : : end subroutine tdep_calc_write_gruneisen
2262 : : !!***
2263 : :
2264 : : !=====================================================================================================
2265 : :
2266 : : !FB!=====================================================================================================
2267 : : !FBsubroutine tdep_calc_lifetime1(Crystal,distance,Eigen2nd,Ifc,Invar,Lattice,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym)
2268 : : !FB
2269 : : !FB type(crystal_t),intent(in) :: Crystal
2270 : : !FB type(Symmetries_type),intent(in) :: Sym
2271 : : !FB type(atdep_dataset_type),intent(in) :: Invar
2272 : : !FB type(Shell_type),intent(in) :: Shell3at
2273 : : !FB type(Lattice_type),intent(in) :: Lattice
2274 : : !FB type(Eigen_type),intent(in) :: Eigen2nd
2275 : : !FB type(Qbz_type),intent(in) :: Qbz
2276 : : !FB type(ifc_type),intent(in) :: Ifc
2277 : : !FB double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
2278 : : !FB double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
2279 : : !FB double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
2280 : : !FB
2281 : : !FB integer :: iq_bz,jq_bz,kq_bz,nmode,ii,tot_count
2282 : : !FB integer :: okp_count,okm_count
2283 : : !FB double precision :: iqbz(3),jqbz(3),kqbz(3),ipjpk(3),imjpk(3)
2284 : : !FB! double complex, allocatable :: Gruneisen(:,:,:)
2285 : : !FB! double complex, allocatable :: Grun_mean(:)
2286 : : !FB integer :: ierr
2287 : : !FB
2288 : : !FB ierr = 0
2289 : : !FB
2290 : : !FB nmode=3*Invar%natom_unitcell
2291 : : !FB tot_count=zero
2292 : : !FB okp_count=zero
2293 : : !FB okm_count=zero
2294 : : !FB do ii=1,3
2295 : : !FB write(ab_out,'(a,3(e15.6,1x),1x)') 'Gmet=',Lattice%gmet(ii,:)
2296 : : !FB write(ab_out,'(a,3(e15.6,1x),1x)') 'Gprim=',Lattice%gprim(ii,:)
2297 : : !FB end do
2298 : : !FB do iq_bz=1,Qbz%nqbz
2299 : : !FB iqbz(:)=Qbz%qbz(:,iq_bz)
2300 : : !FB! iqbz(:)=0.2
2301 : : !FB write(ab_out,'(a,3(e15.6,1x),1x)') 'Qbz1=',iqbz(:)
2302 : : !FB do jq_bz=1,Qbz%nqbz
2303 : : !FB jqbz(:)=Qbz%qbz(:,jq_bz)
2304 : : !FB! jqbz(:)=0.3
2305 : : !FB do kq_bz=1,Qbz%nqbz
2306 : : !FB tot_count=tot_count+1
2307 : : !FB kqbz(:)=Qbz%qbz(:,kq_bz)
2308 : : !FB! kqbz(:)=0.5
2309 : : !FB ipjpk(:)=iqbz(:)+jqbz(:)+kqbz(:)
2310 : : !FB imjpk(:)=iqbz(:)-jqbz(:)+kqbz(:)
2311 : : !FB if (sum(abs(ipjpk(:)-int(ipjpk(:)))).lt.tol4) then
2312 : : !FB write(ab_out,'(a,i4,1x,i4,1x,i4)')'OK for qi+qj+qk=G',iq_bz,jq_bz,kq_bz
2313 : : !FB okp_count=okp_count+1
2314 : : !FB else if (sum(abs(imjpk(:)-int(imjpk(:)))).lt.tol4) then
2315 : : !FB write(ab_out,'(a,i4,1x,i4,1x,i4)')'OK for qi-qj+qk=G',iq_bz,jq_bz,kq_bz
2316 : : !FB okm_count=okm_count+1
2317 : : !FB else
2318 : : !FB cycle
2319 : : !FB end if
2320 : : !FB! do imode=1,nmode
2321 : : !FB! do jmode=1,nmode
2322 : : !FB! do kmode=1,nmode
2323 : : !FB! end do
2324 : : !FB! end do
2325 : : !FB! end do
2326 : : !FB end do
2327 : : !FB end do
2328 : : !FB end do
2329 : : !FB write(ab_out,'(a,i10)') 'Total number of (q1,q2,q3) =',tot_count
2330 : : !FB write(ab_out,'(a,i10)') 'Partial number of (q1+q2+q3) =',okp_count
2331 : : !FB write(ab_out,'(a,i10)') 'Partial number of (q1-q2+q3) =',okm_count
2332 : : !FB
2333 : : !FBend subroutine tdep_calc_lifetime1
2334 : :
2335 : : !=====================================================================================================
2336 : : !FBsubroutine tdep_calc_lifetime2(Crystal,distance,Eigen2nd,Ifc,Invar,Lattice,MPIdata,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym)
2337 : : !FB
2338 : : !FB!Arguments ------------------------------------
2339 : : !FB!scalars
2340 : : !FB integer,intent(in) :: ncid,comm
2341 : : !FB real(dp),intent(in) :: dosdeltae !,dossmear
2342 : : !FB type(gruns_t),intent(in) :: gruns
2343 : : !FB character(len=*),intent(in) :: prefix
2344 : : !FB!arrays
2345 : : !FB
2346 : : !FB type(crystal_t),intent(in) :: Crystal
2347 : : !FB type(Symmetries_type),intent(in) :: Sym
2348 : : !FB type(atdep_dataset_type),intent(in) :: Invar
2349 : : !FB type(Shell_type),intent(in) :: Shell3at
2350 : : !FB type(Lattice_type),intent(in) :: Lattice
2351 : : !FB type(Eigen_type),intent(in) :: Eigen2nd
2352 : : !FB type(MPI_enreg_type), intent(in) :: MPIdata
2353 : : !FB type(Qbz_type),intent(in) :: Qbz
2354 : : !FB type(ifc_type),intent(in) :: Ifc
2355 : : !FB double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
2356 : : !FB double precision,intent(in) :: Phi3_ref(3,3,3,Shell3at%nshell)
2357 : : !FB double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
2358 : : !FB
2359 : : !FB!Local variables-------------------------------
2360 : : !FB!scalars
2361 : : !FB integer,parameter :: master=0,qptopt1=1,bcorr0=0
2362 : : !FB integer :: nprocs,my_rank,iqibz,nqbz,nqibz,ierr,ii,nu,ncerr,nomega,cnt,unt,io
2363 : : !FB real(dp) :: gavg,omega_min,omega_max,v2
2364 : : !FB type(htetra_t) :: tetra
2365 : : !FB character(len=500) :: msg
2366 : : !FB!arrays
2367 : : !FB integer :: qptrlatt(3,3)
2368 : : !FB integer :: ngqpt(3)
2369 : : !FB real(dp) :: shiftq(3,1)
2370 : : !FB real(dp),allocatable :: gvals_qibz(:,:),wvols_qibz(:,:,:),dwdq_qibz(:,:,:)
2371 : : !FB real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:)
2372 : : !FB real(dp),allocatable :: wdt(:,:),wdos(:,:),grdos(:,:),gr2dos(:,:),wibz(:),omega_mesh(:)
2373 : : !FB real(dp),allocatable :: vdos(:,:),v2dos(:,:)
2374 : : !FB real(dp),allocatable :: phdispl_cart_qibz(:,:,:,:,:)
2375 : : !FB
2376 : : !FB! ************************************************************************
2377 : : !FB
2378 : : !FB comm = MPIdata%comm_step
2379 : : !FB nprocs = xmpi_comm_size(comm)
2380 : : !FB my_rank = xmpi_comm_rank(comm)
2381 : : !FB
2382 : : !FB write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Lifetime ',ch10
2383 : : !FB call wrtout(ab_out,msg)
2384 : : !FB
2385 : : !FB ! Generate the q-mesh by finding the IBZ and the corresponding weights.
2386 : : !FB qptrlatt = 0
2387 : : !FB do ii=1,3
2388 : : !FB qptrlatt(ii,ii) = Invar%ngqpt1(:)
2389 : : !FB end do
2390 : : !FB nshiftq=1
2391 : : !FB shiftq(:,:)=0.0
2392 : : !FB
2393 : : !FB ! Get IBZ and BZ.
2394 : : !FB call kpts_ibz_from_kptrlatt(Crystal, qptrlatt, qptopt1, nshiftq, shiftq, &
2395 : : !FB nqibz, qibz, wtq, nqbz, qbz)
2396 : : !FB
2397 : : !FB ! Build tetrahedra
2398 : : !FB tetra = tetra_from_kptrlatt(Crystal, qptopt1, qptrlatt, nshiftq, shiftq, nqibz, qibz, comm, msg, ierr)
2399 : : !FB if (ierr /= 0) ABI_ERROR(msg)
2400 : : !FB
2401 : : !FB ABI_CALLOC(wvols_qibz, (gruns%natom3, gruns%nvols, nqibz))
2402 : : !FB ABI_CALLOC(gvals_qibz, (gruns%natom3, nqibz))
2403 : : !FB ABI_CALLOC(dwdq_qibz, (3, gruns%natom3, nqibz))
2404 : : !FB ABI_CALLOC(phdispl_cart_qibz, (2, gruns%natom3, gruns%natom3, gruns%nvols, nqibz))
2405 : : !FB
2406 : : !FB gavg = zero
2407 : : !FB do iqibz=1,nqibz
2408 : : !FB if (mod(iqibz, nprocs) /= my_rank) cycle ! mpi-parallelism
2409 : : !FB call gruns_fourq(gruns, qibz(:,iqibz), wvols_qibz(:,:,iqibz), gvals_qibz(:,iqibz), &
2410 : : !FB dwdq_qibz(:,:,iqibz), phdispl_cart_qibz(:,:,:,:,iqibz))
2411 : : !FB gavg = gavg + wtq(iqibz) * sum(gvals_qibz(:,iqibz))
2412 : : !FB end do
2413 : : !FB gavg = gavg / gruns%natom3
2414 : : !FB
2415 : : !FB call xmpi_sum(gavg, comm, ierr)
2416 : : !FB call xmpi_sum(wvols_qibz, comm, ierr)
2417 : : !FB call xmpi_sum(gvals_qibz, comm, ierr)
2418 : : !FB call xmpi_sum(dwdq_qibz, comm, ierr)
2419 : : !FB call xmpi_sum(phdispl_cart_qibz, comm, ierr)
2420 : : !FB
2421 : : !FB omega_min = gruns%ifc_vol(gruns%iv0)%omega_minmax(1)
2422 : : !FB omega_max = gruns%ifc_vol(gruns%iv0)%omega_minmax(2)
2423 : : !FB nomega = nint((omega_max - omega_min) / dosdeltae) + 1
2424 : : !FB nomega = max(6, nomega) ! Ensure Simpson integration will be ok
2425 : : !FB
2426 : : !FB ABI_MALLOC(omega_mesh, (nomega))
2427 : : !FB omega_mesh = arth(omega_min, dosdeltae, nomega)
2428 : : !FB omega_max = omega_mesh(nomega)
2429 : : !FB !write(std_out,*)"hello",omega_min,omega_max,dosdeltae,(omega_max-omega_min) / (nomega-1)
2430 : : !FB ABI_MALLOC(wibz, (nqibz))
2431 : : !FB ABI_MALLOC(wdt, (nomega, 2))
2432 : : !FB ABI_CALLOC(wdos, (nomega, 2))
2433 : : !FB ABI_CALLOC(grdos, (nomega, 2))
2434 : : !FB ABI_CALLOC(gr2dos, (nomega, 2))
2435 : : !FB ABI_CALLOC(vdos, (nomega, 2))
2436 : : !FB ABI_CALLOC(v2dos, (nomega, 2))
2437 : : !FB
2438 : : !FB ! Compute DOSes.
2439 : : !FB cnt = 0
2440 : : !FB do iqibz=1,nqibz
2441 : : !FB do nu=1,gruns%natom3
2442 : : !FB cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi-parallelism
2443 : : !FB wibz = wvols_qibz(nu, gruns%iv0, :)
2444 : : !FB call tetra%get_onewk(iqibz,bcorr0,nomega,nqibz,wibz,omega_min,omega_max,one,wdt)
2445 : : !FB wdt = wdt*wtq(iqibz)
2446 : : !FB wdos = wdos + wdt
2447 : : !FB grdos = grdos + wdt * gvals_qibz(nu,iqibz)
2448 : : !FB gr2dos = gr2dos + wdt * gvals_qibz(nu,iqibz) ** 2
2449 : : !FB v2 = sum(dwdq_qibz(1:3,nu,iqibz) ** 2)
2450 : : !FB vdos = vdos + wdt * sqrt(v2)
2451 : : !FB v2dos = v2dos + wdt * v2
2452 : : !FB end do
2453 : : !FB end do
2454 : : !FB
2455 : : !FB call xmpi_sum(wdos, comm, ierr)
2456 : : !FB call xmpi_sum(grdos, comm, ierr)
2457 : : !FB call xmpi_sum(gr2dos, comm, ierr)
2458 : : !FB call xmpi_sum(vdos, comm, ierr)
2459 : : !FB call xmpi_sum(v2dos, comm, ierr)
2460 : : !FB
2461 : : !FB if (my_rank == master) then
2462 : : !FB write(unt,'(a)')'# Phonon density of states, Gruneisen DOS and phonon group velocity DOS'
2463 : : !FB write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
2464 : : !FB write(unt,'(a,i0)')'# Tetrahedron method with nqibz= ',nqibz
2465 : : !FB write(unt,"(a,f8.5)")"# Average Gruneisen parameter:", gavg
2466 : : !FB write(unt,'(5a)') &
2467 : : !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"
2468 : : !FB do io=1,nomega
2469 : : !FB write(unt, "(11es17.8)")omega_mesh(io), &
2470 : : !FB wdos(io,1), grdos(io,1), gr2dos(io,1), vdos(io,1), v2dos(io,1), &
2471 : : !FB wdos(io,2), grdos(io,2), gr2dos(io,2), vdos(io,2), v2dos(io,2)
2472 : : !FB end do
2473 : : !FB end if
2474 : : !FB
2475 : : !FB
2476 : : !FB ABI_FREE(qibz)
2477 : : !FB ABI_FREE(wtq)
2478 : : !FB ABI_FREE(qbz)
2479 : : !FB ABI_FREE(wvols_qibz)
2480 : : !FB ABI_FREE(gvals_qibz)
2481 : : !FB ABI_FREE(dwdq_qibz)
2482 : : !FB ABI_FREE(phdispl_cart_qibz)
2483 : : !FB ABI_FREE(omega_mesh)
2484 : : !FB ABI_FREE(wibz)
2485 : : !FB ABI_FREE(wdt)
2486 : : !FB ABI_FREE(wdos)
2487 : : !FB ABI_FREE(grdos)
2488 : : !FB ABI_FREE(gr2dos)
2489 : : !FB ABI_FREE(v2dos)
2490 : : !FB ABI_FREE(vdos)
2491 : : !FB
2492 : : !FB call tetra%free()
2493 : : !FB
2494 : : !FBend subroutine tdep_calc_lifetime2
2495 : : !!***
2496 : : !====================================================================================================
2497 : :
2498 [ - + - - : 138 : end module m_tdep_model
+ + - + +
- - + - +
- + - + -
+ + - - +
- + - + -
+ - + -
+ ]
2499 : : !!***
|
