Branch data Line data Source code
1 : : !!****m* ABINIT/m_phonons
2 : : !! NAME
3 : : !! m_phonons
4 : : !!
5 : : !! FUNCTION
6 : : !! Module for the phonon density of states.
7 : : !! Container type is defined, and destruction, print subroutines
8 : : !! as well as the central phdos_init
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 1999-2026 ABINIT group (XG, MG, MJV, GMR)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public Licence, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_phonons
26 : :
27 : : use, intrinsic :: iso_c_binding
28 : : use defs_basis
29 : : use m_errors
30 : : use m_xmpi
31 : : use m_abicore
32 : : use m_htetra
33 : : use m_numeric_tools
34 : : use m_cgtools
35 : : use m_crystal
36 : : use m_nctk
37 : : use m_sortph
38 : : use m_ddb
39 : : use netcdf
40 : : use m_supercell
41 : :
42 : :
43 : : use m_fstrings, only : itoa, ftoa, sjoin, ltoa, ktoa, strcat, basename, replace
44 : : use m_matrix, only : matr3inv
45 : : use m_symtk, only : sg_multable
46 : : use m_krank, only : krank_t
47 : : use m_dtset, only : dataset_type
48 : : use m_time, only : cwtime, cwtime_report
49 : : use m_io_tools, only : open_file
50 : : use m_geometry, only : mkrdim, symredcart, phdispl_cart2red
51 : : use m_dynmat, only : gtdyn9, dfpt_phfrq, dfpt_prtph, &
52 : : pheigvec_normalize, massmult_and_breaksym, phdispl_from_eigvec, phangmom_from_eigvec
53 : : use m_atprj, only : atprj_type
54 : : use m_bz_mesh, only : isamek, make_path, kpath_t
55 : : use m_ifc, only : ifc_type
56 : : use m_anaddb_dataset, only : anaddb_dataset_type
57 : : use m_kpts, only : kpts_ibz_from_kptrlatt, get_full_kgrid, kpts_map, kpts_timrev_from_kptopt
58 : : use m_special_funcs, only : bose_einstein
59 : : use m_sort, only : sort_dp
60 : : use m_symfind, only : symanal
61 : :
62 : : implicit none
63 : :
64 : : private
65 : :
66 : : public :: mkphbs ! Compute phonon band structure
67 : : public :: phonons_write_xmgrace ! Write phonons bands in Xmgrace format.
68 : : public :: phonons_write_gnuplot ! Write phonons bands in gnuplot format.
69 : : public :: ifc_mkphbs ! Compute the phonon band structure from the IFC and write data to file(s)
70 : : public :: dfpt_symph ! Determine the symmetry character of the different phonon modes at Gamma
71 : :
72 : : public :: zacharias_supercell_make
73 : : public :: zacharias_supercell_print
74 : : public :: thermal_supercell_make
75 : : public :: thermal_supercell_free
76 : : public :: thermal_supercell_print
77 : : !!***
78 : :
79 : : !!****t* m_phonons/phdos_t
80 : : !! NAME
81 : : !! phdos_t
82 : : !!
83 : : !! FUNCTION
84 : : !! Container for phonon DOS and atom projected contributions
85 : : !!
86 : : !! SOURCE
87 : :
88 : : type,public :: phdos_t
89 : :
90 : : integer :: ntypat
91 : : ! Number of type of atoms.
92 : :
93 : : integer :: natom
94 : : ! Number of atoms is the unit cell.
95 : :
96 : : integer :: nomega
97 : : ! Number of frequency points in DOS mesh.
98 : :
99 : : integer :: nqibz
100 : : ! Number of q-points in the IBZ.
101 : :
102 : : integer :: n_normal_vec_dmm
103 : : ! Number of directions to project the thermal current in the diffuse mismatch model
104 : :
105 : : integer :: prtdos
106 : : ! Option of DOS calculation (1 for Gaussian, 2 for tetrahedrons).
107 : :
108 : : real(dp) :: omega_min
109 : : ! Min frequency for DOS calculation.
110 : :
111 : : real(dp) :: omega_max
112 : : ! Max frequency for DOS calculation.
113 : :
114 : : real(dp) :: omega_step
115 : : ! Frequency step.
116 : :
117 : : real(dp) :: dossmear
118 : : ! Gaussian broadening.
119 : :
120 : : integer :: qptrlatt(3,3) = 0
121 : : ! q-mesh as computed in getkgrid_low
122 : :
123 : : real(dp) :: shiftq(3)
124 : : ! Shigt of Q-mesh computed by getkgrid_low (1 shift is enough)
125 : :
126 : : real(dp),allocatable :: atom_mass(:)
127 : : ! atom_mass(natom)
128 : :
129 : : real(dp),allocatable :: normal_vec_dmm(:,:)
130 : : ! normal_vec_dmm(3,n_normal_vec_dmm)
131 : : ! Cartesian directions to project the thermal current in the diffuse mismatch model
132 : :
133 : : real(dp),allocatable :: omega(:)
134 : : ! omega(nomega)
135 : : ! Frequency grid.
136 : :
137 : : real(dp),allocatable :: phdos(:)
138 : : ! phdos(nomega)
139 : : ! phonon DOS.
140 : :
141 : : real(dp),allocatable :: phdos_int(:)
142 : : ! phdos_int(nomega)
143 : : ! integrated phonon DOS
144 : :
145 : : real(dp),allocatable :: phdos_dmm(:,:)
146 : : ! phdos_dmm(nomega,n_normal_vec_dmm)
147 : : ! phonon DOS including projection of velocity onto a given direction n normal to an interface
148 : : ! in the diffuse mismatch model
149 : :
150 : : real(dp),allocatable :: pjdos(:,:,:)
151 : : ! pjdos(nomega,3,natom)
152 : : ! projected DOS (over atoms and cartesian directions)
153 : :
154 : : real(dp),allocatable :: pjdos_int(:,:,:)
155 : : ! pjdos_int(nomega,3,natom)
156 : : ! Integrated atomic PJDOS along the three cartesian directions.
157 : :
158 : : real(dp),allocatable :: pjdos_type(:,:)
159 : : ! pjdos_type(nomega,ntypat)
160 : : ! phonon DOS contribution arising from a particular atom-type.
161 : :
162 : : real(dp),allocatable :: pjdos_type_int(:,:)
163 : : ! pjdos_type_int(nomega,ntypat)
164 : : ! Integrate phonon DOS contribution arising from a particular atom-type.
165 : :
166 : : real(dp),allocatable :: pjdos_rc_type(:,:,:)
167 : : ! phdos(nomega,3,ntypat)
168 : : ! phonon DOS contribution arising from a particular atom-type
169 : : ! decomposed along the three cartesian directions.
170 : :
171 : : real(dp),allocatable :: msqd_dos_atom(:,:,:,:)
172 : : ! msqd_dos_atom(nomega,3,3,natom)
173 : : ! mean square displacement matrix, frequency dependent like a DOS, tensor in cartesian coords.
174 : : ! allows one to calculate Debye Waller factors by integration with 1/omega
175 : : ! and the Bose Einstein factor
176 : :
177 : : contains
178 : :
179 : : procedure :: print => phdos_print
180 : : procedure :: print_debye => phdos_print_debye
181 : : procedure :: print_msqd => phdos_print_msqd
182 : : procedure :: print_thermo => phdos_print_thermo
183 : : procedure :: free => phdos_free
184 : : procedure :: ncwrite => phdos_ncwrite
185 : : procedure :: init => phdos_init ! Constructor
186 : : end type phdos_t
187 : : !!***
188 : :
189 : : !!****t* m_phonons/phstore_t
190 : : !! NAME
191 : : !! phstore_t
192 : : !!
193 : : !! FUNCTION
194 : : !! This object stores ph eigenvalues and eigenvectors in the IBZ and provides methods
195 : : !! to compute the corresponding quantities in the full BZ using symmetries.
196 : : !! Useful for very intensive loops of q-points in the BZ in which the call to ifc_fourq
197 : : !! may become a significant bottleneck.
198 : : !! Note that IBZ quantities are memory-distributed inside the MPI communicator comm.
199 : : !! and the symmetrization is performed in a non-blocking fashion so that it's possible to overlap
200 : : !! the symmetrization with computations. See sigmaph for usage.
201 : : !!
202 : : !! SOURCE
203 : :
204 : : type,public :: phstore_t
205 : :
206 : : integer :: nqibz
207 : : ! Number of q-points in the IBZ
208 : :
209 : : integer :: comm
210 : : ! MPI communicator used to distribute memory.
211 : :
212 : : integer :: nprocs
213 : : ! Number of MPI procs in comm.
214 : :
215 : : integer :: my_rank
216 : : ! Rank of this MPI proc inside comm
217 : :
218 : : integer :: natom, natom3
219 : :
220 : : logical :: use_ifc_fourq = .False.
221 : : ! Debugging flag. If True, replace symmetrization with call to ifc_fourq.
222 : :
223 : : integer :: requests(2)
224 : : ! MPI requests
225 : :
226 : : integer,allocatable :: qibz_start(:), qibz_stop(:)
227 : : ! (0:%nprocs-1))
228 : : ! Initial and final index of the IBZ qpoint treated by this MPI proc inside comm.
229 : :
230 : : real(dp), contiguous, pointer :: qibz(:,:)
231 : : ! q-points in the IBZ.
232 : :
233 : : real(dp),allocatable :: phfreqs_qibz(:,:)
234 : : ! (natom3, %nqibz))
235 : : ! Ph frequencies in the IBZ
236 : :
237 : : real(dp),allocatable :: pheigvec_qibz(:,:,:,:)
238 : : ! (2, natom3, natom3, %nqibz))
239 : : ! Ph eigenvectors in the IBZ
240 : :
241 : : real(dp),allocatable :: phfrq(:)
242 : : ! (natom3)
243 : : ! Ph frequencies for q in the full BZ
244 : :
245 : : real(dp),allocatable :: displ_cart(:,:,:,:)
246 : : ! displ_cart(2, 3, natom, %natom3)
247 : : ! Ph displacement for q in the full BZ
248 : :
249 : : contains
250 : :
251 : : procedure :: init => phstore_init ! Creation method (allocates memory, initialize data from input vars).
252 : : procedure :: async_rotate => phstore_async_rotate ! Begin non-blocking collective MPI communication to symmetrize stuff
253 : : procedure :: wait => phstore_wait ! Wait from non-blocking MPI BCAST started in phstore_async_rotate,
254 : : ! return ph frequencies and displacements.
255 : : procedure :: free => phstore_free ! Free dynamic memory
256 : : end type phstore_t
257 : : !!***
258 : :
259 : : public :: pheigvec_rotate ! Obtain phonon eigenvectors for q in the BZ from the symmetrical image in the IBZ.
260 : : public :: test_phrotation ! Validate pheigvec_rotate routine.
261 : :
262 : : contains !=====================================================
263 : : !!***
264 : :
265 : : !!****f* m_phonons/phdos_print
266 : : !! NAME
267 : : !! phdos_print
268 : : !!
269 : : !! FUNCTION
270 : : !! Print out phonon DOS (and partial DOS etc) in meV units
271 : : !!
272 : : !! INPUTS
273 : : !! PHdos= container object for phonon DOS
274 : : !! fname=File name for output
275 : : !!
276 : : !! OUTPUT
277 : : !! Only writing.
278 : : !!
279 : : !! SOURCE
280 : :
281 : 79 : subroutine phdos_print(PHdos, fname)
282 : :
283 : : !Arguments ------------------------------------
284 : : class(phdos_t),intent(in) :: PHdos
285 : : character(len=*),intent(in) :: fname
286 : :
287 : : !Local variables-------------------------------
288 : : integer :: io,itype,unt,unt_by_atom,unt_msqd,iatom
289 : : real(dp) :: tens(3,3)
290 : : character(len=500) :: msg, msg_method
291 : : character(len=fnlen) :: fname_by_atom, fname_msqd
292 : : character(len=3) :: unitname
293 : : ! *************************************************************************
294 : :
295 : : ! Use Ha units everywhere
296 : 79 : unitname='Ha'
297 : :
298 : 133 : select case (PHdos%prtdos)
299 : : case (1)
300 : : write(msg_method,'(a,es16.8,2a,i0)')&
301 : 54 : '# Gaussian method with smearing = ',PHdos%dossmear,unitname,', nqibz =',PHdos%nqibz
302 : : case (2)
303 : 25 : write(msg_method,'(a,i0)')'# Tetrahedron method, nqibz= ',PHdos%nqibz
304 : : case default
305 [ + + - ]: 79 : ABI_ERROR(sjoin(" Wrong prtdos: ",itoa(PHdos%prtdos)))
306 : : end select
307 : :
308 : : ! Open external file and write results
309 [ - + ]: 79 : if (open_file(fname,msg,newunit=unt,form="formatted",action="write") /= 0) then
310 : 0 : ABI_ERROR(msg)
311 : : end if
312 : 79 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom type projected DOS'
313 : 79 : call wrtout(unt,msg)
314 : 79 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
315 : 79 : call wrtout(unt,msg)
316 : 79 : call wrtout(unt,msg_method)
317 : 79 : write(msg,'(5a)')'# ',ch10,'# omega PHDOS INT_PHDOS PJDOS[atom_type=1] INT_PJDOS[atom_type=1] ... ',ch10,'# '
318 : 79 : call wrtout(unt,msg)
319 [ + + ]: 154639 : do io=1,PHdos%nomega
320 : 154560 : write(unt,'(3es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io),PHdos%phdos_int(io)
321 [ + + ]: 368390 : do itype=1,PHdos%ntypat
322 : 368390 : write(unt,'(2es17.8,2x)',advance='NO')PHdos%pjdos_type(io,itype),PHdos%pjdos_type_int(io,itype)
323 : : end do
324 : 154639 : write(unt,*)
325 : : end do
326 : 79 : close(unt)
327 : :
328 [ + - ]: 79 : fname_by_atom = trim(fname) // "_by_atom"
329 [ - + ]: 79 : if (open_file(fname_by_atom,msg,newunit=unt_by_atom,form="formatted",action="write") /= 0) then
330 : 0 : ABI_ERROR(msg)
331 : : end if
332 : 79 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom projected DOS'
333 : 79 : call wrtout(unt_by_atom,msg)
334 : 79 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
335 : 79 : call wrtout(unt_by_atom,msg)
336 : 79 : call wrtout(unt_by_atom,msg_method)
337 : 79 : write(msg,'(5a)')'# ',ch10,'# omega PHDOS PJDOS[atom=1] PJDOS[atom=2] ... ',ch10,'# '
338 : 79 : call wrtout(unt_by_atom,msg)
339 [ + + ]: 154639 : do io=1,PHdos%nomega
340 : 154560 : write(unt_by_atom,'(2es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io)
341 [ + + ]: 531829 : do iatom=1,PHdos%natom
342 [ + + ]: 1663636 : write(unt_by_atom,'(1es17.8,2x)',advance='NO') sum(PHdos%pjdos(io,1:3,iatom))
343 : : end do
344 : 154639 : write(unt_by_atom,*)
345 : : end do
346 : 79 : close(unt_by_atom)
347 : :
348 [ + - ]: 79 : fname_msqd = trim(fname) // "_msqd"
349 [ - + ]: 79 : if (open_file(fname_msqd,msg,newunit=unt_msqd,form="formatted",action="write") /= 0) then
350 : 0 : ABI_ERROR(msg)
351 : : end if
352 : 79 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states weighted msq displacement matrix (set to zero below 1e-12)'
353 : 79 : call wrtout(unt_msqd,msg)
354 : 79 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in bohr^2 states/',unitname
355 : 79 : call wrtout(unt_msqd,msg)
356 : 79 : call wrtout(unt_msqd,msg_method)
357 : 79 : write(msg,'(5a)')'# ',ch10,'# omega MSQDisp[atom=1, xx, yy, zz, yz, xz, xy] MSQDisp[atom=2, xx, yy,...] ... ',ch10,'# '
358 : 79 : call wrtout(unt_msqd,msg)
359 [ + + ]: 154639 : do io=1,PHdos%nomega
360 : 154560 : write(unt_msqd,'(2es17.8)',advance='NO')PHdos%omega(io)
361 [ + + ]: 531829 : do iatom=1,PHdos%natom
362 [ + + + + ]: 4904497 : tens = PHdos%msqd_dos_atom(io,:,:,iatom)
363 [ + + + + : 4904497 : where (abs(tens) < tol12)
+ + ]
364 : : tens = zero
365 : : end where
366 : : write(unt_msqd,'(6es17.8,2x)',advance='NO') &
367 : 377269 : tens(1,1), tens(2,2), tens(3,3), &
368 : 909098 : tens(2,3), tens(1,3), tens(1,2)
369 : : end do
370 : 154639 : write(unt_msqd,*)
371 : : end do
372 : 79 : close(unt_msqd)
373 : :
374 : 79 : end subroutine phdos_print
375 : : !!***
376 : :
377 : : !----------------------------------------------------------------------
378 : :
379 : : !****f* m_phonons/phdos_print_debye
380 : : !!
381 : : !! NAME
382 : : !! phdos_print_debye
383 : : !!
384 : : !! FUNCTION
385 : : !! Print out global Debye temperature, force constant, etc... from phonon DOS
386 : : !!
387 : : !! INPUTS
388 : : !! phonon_dos= container object for phonon DOS
389 : : !! ucvol = unit cell volume
390 : : !!
391 : : !! OUTPUT
392 : : !! Only writing.
393 : : !!
394 : : !! SOURCE
395 : :
396 : 11 : subroutine phdos_print_debye(PHdos, ucvol)
397 : :
398 : : !Arguments ------------------------------------
399 : : class(phdos_t),intent(in) :: PHdos
400 : : real(dp), intent(in) :: ucvol
401 : :
402 : : !Local variables-------------------------------
403 : : integer :: io, iomax, iomin
404 : : real(dp) :: avgom2dos, avgspeedofsound, debyefreq, meanfreq, meanfreq2
405 : : character(len=500) :: msg
406 : : !arrays
407 : : integer :: units(2)
408 : 11 : real(dp), allocatable :: om2dos(:), om1dos(:), intdos(:)
409 : : ! *************************************************************************
410 : :
411 [ + + ]: 33 : units = [std_out, ab_out]
412 : :
413 : : ! average speed of sound: coefficient of omega^2 in the DOS is = Volume / 2 pi^2 hbar^3 v_s^3
414 : : ! first find how far out we can fit with a parabola
415 [ - + + - : 40817 : ABI_CALLOC(om2dos, (PHdos%nomega))
+ - + - -
+ - + +
+ ]
416 [ - + + - : 40806 : ABI_CALLOC(om1dos, (PHdos%nomega))
- + - + +
+ ]
417 [ - + + - : 40806 : ABI_CALLOC(intdos, (PHdos%nomega))
- + - + +
+ ]
418 : 40795 : avgom2dos = zero
419 [ + + ]: 40795 : do io=1,PHdos%nomega
420 [ + + ]: 40795 : if (abs(PHdos%omega(io)) > 1.e-8) then
421 : 40780 : om1dos(io) = PHdos%phdos(io) / PHdos%omega(io)
422 : 40780 : om2dos(io) = PHdos%phdos(io) / PHdos%omega(io)**2
423 : : end if
424 : : end do
425 : :
426 : : ! integrate dos / omega
427 : 11 : call simpson_int(PHdos%nomega,PHdos%omega_step,om1dos,intdos)
428 : 11 : meanfreq = intdos(PHdos%nomega)
429 : :
430 : : ! integrate dos / omega^2
431 [ + + ]: 40795 : intdos = zero
432 : 11 : call simpson_int(PHdos%nomega,PHdos%omega_step,om2dos,intdos)
433 : 11 : meanfreq2 = intdos(PHdos%nomega)
434 : :
435 : 11 : iomin = 1; iomax = PHdos%nomega
436 [ + - ]: 10345 : do io = 1, PHdos%nomega
437 : : ! skip eventual negative frequency modes
438 [ + + ]: 10345 : if (PHdos%omega(io) <= tol10) then
439 : : iomin = io
440 : : cycle
441 : : end if
442 : :
443 : : ! accumulate dos * om^2 to make an average
444 : 11 : avgom2dos = avgom2dos + om2dos(io)
445 : : ! first deviation from initial value of more than 10 percent
446 [ - + ]: 11 : if (abs(one-om2dos(iomin)/om2dos(io)) > 0.1_dp) then
447 : : iomax = io; exit
448 : : end if
449 : : end do
450 : :
451 : 11 : avgom2dos = avgom2dos / (iomax-iomin+1)
452 : : ! this value is also useful for partial atomic DOS, related to kinetic energy and Force constant in Moessbauer
453 : :
454 : 11 : avgspeedofsound = (ucvol / 2 / pi**2 / avgom2dos)**third
455 : 11 : write (msg,'(a,E20.10,3a,F16.4,2a)') ' Average speed of sound: ', avgspeedofsound, ' (at units) ',ch10,&
456 : 22 : '- = ', avgspeedofsound * Bohr_Ang * 1.d-13 / Time_Sec, ' [km/s]',ch10
457 : 11 : call wrtout(units, msg)
458 : :
459 : : ! Debye frequency = vs * (6 pi^2 natom / ucvol)**1/3
460 : 11 : debyefreq = avgspeedofsound * (six*pi**2/ucvol)**(1./3.)
461 : 11 : write (msg,'(a,E20.10,3a,E20.10,a)') ' Debye frequency from DOS: ', debyefreq, ' (Ha) ',ch10,&
462 : 22 : '- = ', debyefreq*Ha_THz, ' (THz)'
463 : 11 : call wrtout(units, msg)
464 : :
465 : : ! Debye temperature = hbar * Debye frequency / kb
466 : 11 : write (msg,'(a,E20.10,2a)') '-Debye temperature from DOS: ', debyefreq*Ha_K, ' (K)', ch10
467 : 11 : call wrtout(units, msg)
468 : :
469 [ - + ]: 11 : ABI_FREE(om2dos)
470 [ - + ]: 11 : ABI_FREE(om1dos)
471 : 11 : ABI_FREE(intdos)
472 : :
473 : 11 : end subroutine phdos_print_debye
474 : : !!***
475 : :
476 : : !----------------------------------------------------------------------
477 : :
478 : : !****f* m_phonons/phdos_print_thermo
479 : : !! NAME
480 : : !! phdos_print_thermo
481 : : !!
482 : : !! FUNCTION
483 : : !! Print out global thermodynamic quantities based on DOS
484 : : !! Only master node should call this routine.
485 : : !!
486 : : !! INPUTS
487 : : !! phonon_dos= container object for phonon DOS
488 : : !! ucvol = unit cell volume
489 : : !!
490 : : !! OUTPUT
491 : : !! Only writing.
492 : : !!
493 : : !! SOURCE
494 : :
495 : 11 : subroutine phdos_print_thermo(PHdos, fname, ntemper, tempermin, temperinc)
496 : :
497 : : !Arguments ------------------------------------
498 : : class(phdos_t),intent(in) :: PHdos
499 : : integer, intent(in) :: ntemper
500 : : real(dp), intent(in) :: tempermin, temperinc
501 : : character(len=*),intent(in) :: fname
502 : :
503 : : !Local variables-------------------------------
504 : : integer :: iomega, itemper, tunt
505 : : character(len=500) :: msg
506 : : real(dp) :: wover2t, ln2shx, cothx, invsinh2
507 : : real(dp) :: tmp, domega
508 : : !arrays
509 : 11 : real(dp), allocatable :: free(:), energy(:), entropy(:), spheat(:),wme(:)
510 : :
511 : : ! *********************************************************************
512 : :
513 : : ! Allocate and put zeroes for F, E, S, Cv
514 [ - + + - : 143 : ABI_CALLOC(free, (ntemper))
+ - + - -
+ - + +
+ ]
515 [ - + + - : 132 : ABI_CALLOC(energy, (ntemper))
- + - + +
+ ]
516 [ - + + - : 132 : ABI_CALLOC(entropy, (ntemper))
- + - + +
+ ]
517 [ - + + - : 132 : ABI_CALLOC(spheat, (ntemper))
- + - + +
+ ]
518 [ - + + - : 132 : ABI_CALLOC(wme, (ntemper))
- + - + +
+ ]
519 : :
520 : : ! open THERMO file
521 [ - + ]: 11 : if (open_file(fname, msg, newunit=tunt, form="formatted", action="write") /= 0) then
522 : 0 : ABI_ERROR(msg)
523 : : end if
524 : :
525 : : write(msg, '(3a)' )&
526 : 11 : ' phdos_print_thermo: thermodynamic functions calculated from prtdos DOS (not histogram)',ch10,&
527 : 22 : ' see THERMO output file ...'
528 : 11 : call wrtout(std_out,msg)
529 : :
530 : : ! print header
531 : 11 : write(tunt,'(a,a)') ch10,&
532 : 22 : ' # At T F(J/mol-c) E(J/mol-c) S(J/(mol-c.K)) C(J/(mol-c.K)) Omega_mean(cm-1) from prtdos DOS'
533 : 11 : write(tunt, "(a)")' # (A mol-c is the abbreviation of a mole-cell, that is, the'
534 : 11 : write(tunt, "(a)")' # number of Avogadro times the atoms in a unit cell)'
535 : :
536 : 11 : domega = phdos%omega_step
537 : :
538 [ + + ]: 121 : do itemper=1,ntemper
539 : : ! The temperature (tmp) is given in Ha
540 : 110 : tmp=(tempermin+temperinc*dble(itemper-1))*kb_HaK
541 : :
542 [ + + ]: 407950 : do iomega=1,PHdos%nomega
543 [ + + ]: 407840 : if (abs(PHdos%phdos(iomega)) < 1.e-200_dp) cycle
544 : :
545 : : ! wover2t= hbar*w / 2kT dimensionless
546 [ + - ]: 375470 : wover2t = zero; if(tmp > tol14) wover2t=PHdos%omega(iomega)*half/tmp
547 : : ! should not be much of a problem for the log, but still put a check.
548 [ + + + - ]: 375470 : ln2shx=zero; if (wover2t > tol16 .and. wover2t < 100.0_dp) ln2shx=log(two * sinh(wover2t))
549 [ + - ]: 281080 : cothx=zero; if (wover2t > tol16) cothx=one/tanh(wover2t)
550 [ + - ]: 281080 : invsinh2=zero; if (wover2t > tol16 .and. wover2t < 100.0_dp) invsinh2=one/sinh(wover2t)**2
551 : :
552 : : ! This matches the equations published in Lee & Gonze, PRB 51, 8610 (1995) [[cite:Lee1995]]
553 : 375470 : free(itemper) = free(itemper) + PHdos%phdos(iomega)*tmp*ln2shx
554 : 375470 : energy(itemper) = energy(itemper) + PHdos%phdos(iomega)*half*PHdos%omega(iomega)*cothx
555 : 375470 : spheat(itemper) = spheat(itemper) + PHdos%phdos(iomega)*wover2t**2 * invsinh2
556 : 375470 : entropy(itemper)= entropy(itemper) + PHdos%phdos(iomega)*(wover2t*cothx - ln2shx)
557 : 407950 : wme(itemper) = wme(itemper) + PHdos%phdos(iomega)*PHdos%omega(iomega)*wover2t**2 * invsinh2
558 : : end do ! iomega
559 : :
560 : : ! suppose homogeneous omega grid and multiply by domega
561 : 110 : free(itemper) = free(itemper) * domega
562 : 110 : energy(itemper) = energy(itemper) * domega
563 : 110 : entropy(itemper)= entropy(itemper) * domega
564 : 110 : spheat(itemper) = spheat(itemper) * domega
565 : 110 : wme(itemper) = wme(itemper) * domega
566 : :
567 [ + - ]: 110 : if (abs(spheat(itemper))>tol8) wme(itemper)=wme(itemper)/spheat(itemper)
568 : :
569 : : ! do the printing to file
570 : 110 : write(tunt,'(es11.3,5es15.7)') tmp/kb_HaK,&
571 : 110 : Ha_J*Avogadro*free(itemper),&
572 : 110 : Ha_J*Avogadro*energy(itemper),&
573 : 110 : Ha_J*Avogadro*kb_HaK*entropy(itemper),&
574 : 110 : Ha_J*Avogadro*kb_HaK*spheat(itemper),&
575 : 231 : wme(itemper)*Ha_cmm1
576 : : end do ! itemper
577 : :
578 : 11 : close(tunt)
579 : :
580 [ - + ]: 11 : ABI_FREE(free)
581 [ - + ]: 11 : ABI_FREE(energy)
582 [ - + ]: 11 : ABI_FREE(entropy)
583 [ - + ]: 11 : ABI_FREE(spheat)
584 [ - + ]: 11 : ABI_FREE(wme)
585 : :
586 : 11 : end subroutine phdos_print_thermo
587 : : !!***
588 : : !----------------------------------------------------------------------
589 : :
590 : : !!****f* m_phonons/phdos_free
591 : : !! NAME
592 : : !! phdos_free
593 : : !!
594 : : !! FUNCTION
595 : : !! Free memory
596 : : !!
597 : : !! SOURCE
598 : :
599 : 177 : subroutine phdos_free(PHdos)
600 : :
601 : : !Arguments -------------------------------
602 : : class(phdos_t),intent(inout) ::PHdos
603 : : ! *************************************************************************
604 : :
605 [ + - ]: 177 : ABI_SFREE(PHdos%atom_mass)
606 [ + - ]: 177 : ABI_SFREE(PHdos%normal_vec_dmm)
607 [ + - ]: 177 : ABI_SFREE(PHdos%omega)
608 [ + - ]: 177 : ABI_SFREE(PHdos%phdos)
609 [ + - ]: 177 : ABI_SFREE(PHdos%phdos_int)
610 [ + - ]: 177 : ABI_SFREE(PHdos%phdos_dmm)
611 [ + - ]: 177 : ABI_SFREE(PHdos%pjdos)
612 [ + - ]: 177 : ABI_SFREE(PHdos%pjdos_int)
613 [ + - ]: 177 : ABI_SFREE(PHdos%pjdos_type)
614 [ + - ]: 177 : ABI_SFREE(PHdos%pjdos_type_int)
615 [ + - ]: 177 : ABI_SFREE(PHdos%pjdos_rc_type)
616 [ + - ]: 177 : ABI_SFREE(PHdos%msqd_dos_atom)
617 : :
618 : 177 : end subroutine phdos_free
619 : : !!***
620 : :
621 : : !--------------------------------------------------------------------------
622 : :
623 : : !!****f* m_phonons/phdos_malloc
624 : : !! NAME
625 : : !! phdos_malloc
626 : : !!
627 : : !! FUNCTION
628 : : !! Allocate memory
629 : : !!
630 : : !! INPUTS
631 : : !!
632 : : !! OUTPUT
633 : : !!
634 : : !! SOURCE
635 : :
636 [ + - ]: 177 : subroutine phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos, n_normal_vec_dmm)
637 : :
638 : : ! Arguments ------------------------------------------------------
639 : : class(phdos_t),intent(out) :: phdos
640 : : type(crystal_t),intent(in) :: crystal
641 : : type(ifc_type),intent(in) :: ifc
642 : : integer,intent(in) :: prtdos
643 : : integer,intent(in),optional :: n_normal_vec_dmm
644 : : real(dp),intent(in) :: dosdeltae,dossmear
645 : : real(dp),intent(in) :: wminmax(2)
646 : :
647 : : !Local variables -------------------------
648 : : integer :: io
649 : : ! *********************************************************************
650 : :
651 : 177 : phdos%ntypat = crystal%ntypat
652 : 177 : phdos%natom = crystal%natom
653 : 177 : phdos%prtdos = prtdos
654 : 177 : phdos%dossmear = dossmear
655 : 177 : phdos%omega_step = dosdeltae
656 : 177 : phdos%n_normal_vec_dmm = 7
657 [ - + ]: 177 : if (present(n_normal_vec_dmm)) then
658 : 0 : phdos%n_normal_vec_dmm = n_normal_vec_dmm
659 : : end if
660 : : ! Use values stored in ifc (obtained with ab-initio q-mesh + pad)
661 [ + + ]: 177 : if (wminmax(2) > wminmax(1)) then
662 : 68 : phdos%omega_min = wminmax(1)
663 : 68 : phdos%omega_max = wminmax(2)
664 : : else
665 : 109 : phdos%omega_min = ifc%omega_minmax(1)
666 : 109 : phdos%omega_max = ifc%omega_minmax(2)
667 : : end if
668 : : ! Must be consistent with mesh computed in tetra routines!
669 : 177 : phdos%nomega = nint((phdos%omega_max - phdos%omega_min) / phdos%omega_step) + 1
670 : : ! Ensure Simpson integration will be ok
671 : 177 : phdos%nomega = max(6, phdos%nomega)
672 : :
673 : : ! Build frequency mesh.
674 [ - + + - : 531 : ABI_MALLOC(phdos%omega, (phdos%nomega))
+ - - + -
+ - + ]
675 [ + + ]: 205541 : do io=1,phdos%nomega
676 : 205541 : phdos%omega(io) = phdos%omega_min + phdos%omega_step * (io - 1)
677 : : end do
678 : 177 : phdos%omega_min = phdos%omega(1)
679 : 177 : phdos%omega_max = phdos%omega(phdos%nomega)
680 : :
681 : : ! Allocate arrays that depend on nomega and set them to zero.
682 [ - + + - : 5487 : ABI_CALLOC(phdos%normal_vec_dmm, (3, phdos%n_normal_vec_dmm))
+ - + - -
+ - + - +
+ + + + ]
683 [ - + + - : 205718 : ABI_CALLOC(phdos%phdos, (phdos%nomega))
- + - + -
+ + + ]
684 [ - + + - : 205718 : ABI_CALLOC(phdos%phdos_int, (phdos%nomega))
- + - + -
+ + + ]
685 [ - + + - : 1439495 : ABI_CALLOC(phdos%phdos_dmm, (phdos%nomega,phdos%n_normal_vec_dmm))
- + + - +
- + - - +
- + - + +
+ + + ]
686 [ - + + - : 1442517 : ABI_CALLOC(phdos%pjdos, (phdos%nomega, 3, crystal%natom))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
687 [ - + + - : 1442340 : ABI_CALLOC(phdos%pjdos_int, (phdos%nomega, 3, crystal%natom))
- + + - -
+ - + - +
+ + + + +
+ ]
688 [ - + + - : 4326690 : ABI_CALLOC(phdos%msqd_dos_atom, (phdos%nomega, 3, 3, crystal%natom))
+ - - + +
- + - + -
- + - + -
+ + + + +
+ + + + ]
689 [ - + + - : 531 : ABI_MALLOC(phdos%atom_mass, (crystal%natom))
+ - + - -
+ - + -
+ ]
690 [ + - + + : 1110 : phdos%atom_mass = crystal%amu(crystal%typat(:)) * amu_emass
+ - - + -
- + + ]
691 : :
692 : 177 : end subroutine phdos_malloc
693 : : !!***
694 : :
695 : : !---------------------------------------------------------------
696 : :
697 : : !!****f* m_phonons/phdos_init
698 : : !! NAME
699 : : !! phdos_init
700 : : !!
701 : : !! FUNCTION
702 : : !! Calculate the phonon density of states as well as
703 : : !! the contributions associated to the different types of atoms in the unit cell.
704 : : !! Two methods are implemented: gaussian method and linear interpolation based on tetrahedra.
705 : : !!
706 : : !! INPUTS
707 : : !! ifc<ifc_type>=Interatomic force constants
708 : : !! crystal<crystal_t>=Info on the crystalline structure.
709 : : !! prtdos=1 for gaussian method, 2 for tetrahedra.
710 : : !! dosdeltae=Step of frequency mesh.
711 : : !! dossmear=Gaussian broadening, used if prtdos==1.
712 : : !! dos_ngqpt(3)=Divisions of the q-mesh used for computing the DOS
713 : : !! nqshift=Number of shifts in Q-mesh
714 : : !! dos_qshift(3, nqshift)=Shift of the q-mesh.
715 : : !! prefix=Prefix for PHBIZ output file. Empty string to deactivate output.
716 : : !! comm=MPI communicator.
717 : : !! prtout=write info to the output and log files.
718 : : !! [exclude_gamma_acoustic]=if .true., the three acoustic modes at Gamma (rigid translations, at zero
719 : : !! frequency) are not included in the DOS (gaussian method only). With a q-mesh containing Gamma,
720 : : !! their smeared contribution at w~0 otherwise makes the integrals of the DOS weighted by 1/w or 1/w**2
721 : : !! (mean square displacement) diverge. Default: .false.
722 : : !!
723 : : !! OUTPUT
724 : : !! phdos<phdos_t>=Container with phonon DOS, IDOS and atom-projected DOS.
725 : : !! count_wminmax(2)=Number of (interpolated) phonon frequencies that are outside
726 : : !! input range (see wminmax). Client code can use count_wminmax and wminmax to
727 : : !! enlarge the mesh and call the routine again to recompute the DOS
728 : : !! if all frequencies should be included.
729 : : !!
730 : : !! SIDE EFFECTS
731 : : !! wminmax(2)=
732 : : !! In input: min and max value of frequency mesh. Used only if minmax(2) > minmax(1)
733 : : !! else values are taken from ifc%omega_minmax (computed from ab-initio mesh + pad)
734 : : !! In output: min and max frequency obtained after interpolating the IFCs on the dense q-mesh dos_ngqpt
735 : : !!
736 : : !! SOURCE
737 : :
738 : 177 : subroutine phdos_init(phdos, crystal, ifc, prtdos, dosdeltae_in, dossmear, dos_ngqpt, nqshft, dos_qshift, prefix, &
739 : : wminmax, count_wminmax, comm, dos_maxmode, prtout, exclude_gamma_acoustic)
740 : :
741 : : !Arguments -------------------------------
742 : : !scalars
743 : : class(phdos_t),intent(out) :: phdos
744 : : integer,intent(in) :: prtdos,nqshft,comm
745 : : real(dp),intent(in) :: dosdeltae_in,dossmear
746 : : character(len=*),intent(in) :: prefix
747 : : type(crystal_t),intent(in) :: crystal
748 : : type(ifc_type),intent(in) :: ifc
749 : : integer, optional, intent(in) :: dos_maxmode
750 : : logical, optional, intent(in) :: exclude_gamma_acoustic
751 : : logical, optional, intent(in) :: prtout
752 : : !arrays
753 : : integer,intent(in) :: dos_ngqpt(3)
754 : : integer,intent(out) :: count_wminmax(2)
755 : : real(dp),intent(in) :: dos_qshift(3,nqshft)
756 : : real(dp),intent(inout) :: wminmax(2)
757 : :
758 : : !Local variables -------------------------
759 : : !scalars
760 : : integer,parameter :: bcorr0 = 0, master = 0
761 : : integer :: iat,jat,idir,imode,io,iq_ibz,itype, my_qptopt, my_nsym
762 : : integer :: nqbz,ierr,natom,nomega,jdir, isym, nprocs, my_rank, ncid, icart
763 : : logical :: refine_dosdeltae
764 : : real(dp),parameter :: max_occ1=one, gaussmaxarg = sqrt(-log(1.d-90)), max_smallq = 0.0625_dp
765 : : real(dp) :: nsmallq,gaussfactor,gaussprefactor,normq,debyefreq,rtmp
766 : : real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all
767 : : real(dp) :: dosdeltae, phdos_int, projfact
768 : : character(len=500) :: msg
769 : : character(len=80) :: errstr
770 : 177 : type(htetra_t) :: htetra_q
771 : : !arrays
772 : : integer :: in_qptrlatt(3,3),new_qptrlatt(3,3), units(2)
773 : : integer :: dos_maxmode_, imin, ii
774 : : logical :: prtout_
775 : 177 : integer,allocatable :: bz2ibz_smap(:,:), bz2ibz(:)
776 : : real(dp) :: veloc(3), speedofsound(3),speedofsound_(3)
777 : 354 : real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom)
778 : 531 : real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom),phfrq(3*Crystal%natom),phangmom(3,3*Crystal%natom)
779 : 354 : real(dp) :: abs_phfrq(3*Crystal%natom)
780 : 354 : logical :: skip_mode(3*Crystal%natom)
781 : 354 : real(dp) :: veloc_1q(3,3*Crystal%natom)
782 : : real(dp) :: qlatt(3,3),rlatt(3,3), msqd_atom_tmp(3,3),temp_33(3,3)
783 : 354 : real(dp) :: symcart(3,3,crystal%nsym), syme2_xyza(3, crystal%natom)
784 : 177 : real(dp),allocatable :: full_eigvec(:,:,:,:,:),full_phfrq(:,:),full_phangmom(:,:,:),new_shiftq(:,:), full_veloc(:,:,:)
785 : 177 : real(dp),allocatable :: qbz(:,:),qibz(:,:),tmp_phfrq(:) !, work_msqd(:,:,:,:)
786 [ + - ]: 177 : real(dp),allocatable :: wtq_ibz(:),xvals(:), gvals_wtq(:), wdt(:,:), energies(:)
787 : : ! *********************************************************************
788 : :
789 : : DBG_ENTER("COLL")
790 : :
791 : 177 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
792 [ + + ]: 531 : units = [std_out, ab_out]
793 : 177 : prtout_ = .true.
794 [ + + ]: 177 : if (present(prtout)) then
795 : 92 : prtout_ = prtout
796 : : end if
797 : :
798 : : ! Consistency check.
799 [ - + ]: 177 : if (all(prtdos /= [1, 2])) then
800 : 0 : ABI_BUG(sjoin('prtdos should be 1 or 2, but received', itoa(prtdos)))
801 : : end if
802 : 177 : dosdeltae = dosdeltae_in; refine_dosdeltae = .false.
803 [ - + ]: 177 : if (dosdeltae <= zero) then
804 : 0 : dosdeltae = -dosdeltae; refine_dosdeltae = .true.
805 [ # # ]: 0 : ABI_CHECK(nprocs == 1, "refine_dosdeltae cannot be used with nprocs > 1")
806 : : end if
807 [ + + - + ]: 177 : if (prtdos == 1 .and. dossmear <= zero) then
808 : 0 : ABI_BUG(sjoin('dossmear should be positive but received', ftoa(dossmear)))
809 : : end if
810 : :
811 : 177 : call cwtime(cpu_all, wall_all, gflops_all, "start")
812 : :
813 : : ! Get symmetries in cartesian coordinates
814 [ + + ]: 7315 : do isym=1,crystal%nsym
815 : 7315 : call symredcart(crystal%rprimd,crystal%gprimd,symcart(:,:,isym),crystal%symrel(:,:,isym))
816 : : end do
817 : :
818 : 177 : natom = crystal%natom
819 : 177 : call phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos)
820 : 177 : nomega = phdos%nomega
821 : :
822 [ + + + + ]: 4713 : veloc_1q = zero
823 : : ! TODO: add input variable and adapt these directions to be read in from input file
824 [ + - ]: 177 : if (phdos%n_normal_vec_dmm >= 7) then
825 [ + + ]: 708 : phdos%normal_vec_dmm(:,1) = [one, zero, zero]
826 [ + + ]: 708 : phdos%normal_vec_dmm(:,2) = [zero, one, zero]
827 [ + + ]: 708 : phdos%normal_vec_dmm(:,3) = [zero, zero, one]
828 [ + + ]: 708 : phdos%normal_vec_dmm(:,4) = [one, one, zero]
829 [ + + ]: 708 : phdos%normal_vec_dmm(:,5) = [one, zero, one]
830 [ + + ]: 708 : phdos%normal_vec_dmm(:,6) = [zero, one, one]
831 [ + + ]: 708 : phdos%normal_vec_dmm(:,7) = [one, one, one]
832 : : end if
833 [ - + + - : 531 : ABI_MALLOC(gvals_wtq, (nomega))
+ - + - -
+ - + ]
834 [ - + + - : 354 : ABI_MALLOC(xvals, (nomega))
- + - + ]
835 : :
836 : : ! Parameters defining the gaussian approximant.
837 [ + + ]: 177 : if (prtdos == 1) then
838 : : ! TODO: use gaussian and update reference files.
839 : 101 : gaussprefactor = one / (dossmear * sqrt(two_pi))
840 : 101 : gaussfactor = one / (sqrt2 * dossmear)
841 : 101 : write(msg, '(4a,f8.5,2a,f8.5,a,i0)') ch10, &
842 : 101 : ' phdos_init: calculating phonon DOS using gaussian method:', ch10, &
843 : 101 : ' gaussian smearing [meV] = ', dossmear * Ha_meV, ch10, &
844 : 202 : ' frequency step [meV] = ', phdos%omega_step * Ha_meV, ", nomega = ",phdos%nomega
845 [ - + ]: 76 : else if (prtdos == 2) then
846 : 76 : write(msg, '(4a,f8.5,a,i0)') ch10, &
847 : 76 : ' phdos_init: calculating phonon DOS using tetrahedron method:', ch10, &
848 : 152 : ' frequency step [meV] = ',phdos%omega_step * Ha_meV, ", nomega = ",phdos%nomega
849 : : end if
850 [ + + ]: 177 : if (prtout_) call wrtout(std_out, msg)
851 : :
852 : : ! This call will set %nqibz and IBZ and BZ arrays
853 : 177 : in_qptrlatt = 0; in_qptrlatt(1, 1) = dos_ngqpt(1); in_qptrlatt(2, 2) = dos_ngqpt(2); in_qptrlatt(3, 3) = dos_ngqpt(3)
854 : :
855 : 177 : my_qptopt = 1
856 : : !my_qptopt = 3 ! This to deactivate the use of symmetries for debugging purposes.
857 : : call kpts_ibz_from_kptrlatt(crystal, in_qptrlatt, my_qptopt, nqshft, dos_qshift, &
858 : 177 : phdos%nqibz, qibz, wtq_ibz, nqbz, qbz, new_kptrlatt=new_qptrlatt, new_shiftk=new_shiftq, bz2ibz=bz2ibz_smap)
859 : :
860 : : my_nsym = crystal%nsym; if (my_qptopt == 3) my_nsym = 1
861 : :
862 [ + + + + ]: 2301 : phdos%qptrlatt = new_qptrlatt
863 [ + + ]: 708 : phdos%shiftq(:) = new_shiftq(:, 1) ! only one shift in output
864 : :
865 [ + - ]: 177 : if (my_rank == master) then
866 [ + - ]: 177 : write(msg, "(3a, i0)")" DOS ngqpt: ", trim(ltoa(dos_ngqpt)), ", qptopt: ", my_qptopt
867 [ + + ]: 177 : if (prtout_) call wrtout(std_out, msg)
868 : 177 : write(msg, "(2(a, i0))")" Number of q-points in the IBZ: ", phdos%nqibz, ", number of MPI processes: ", nprocs
869 [ + + ]: 177 : if (prtout_) call wrtout(std_out, msg)
870 : : end if
871 : : !call cwtime_report(" kpts_ibz_from_kptrlatt", cpu, wall, gflops)
872 : :
873 [ + + ]: 177 : if (prtdos == 2) then
874 : : ! Prepare tetrahedron method including workspace arrays.
875 : : ! Convert kptrlatt to double and invert, qlatt here refer to the shortest qpt vectors
876 [ + + + + ]: 988 : rlatt = new_qptrlatt; call matr3inv(rlatt, qlatt)
877 : :
878 [ - + + - : 228 : ABI_MALLOC(bz2ibz, (nqbz))
+ - + - -
+ - + ]
879 [ + - - + : 472734 : bz2ibz = bz2ibz_smap(1,:)
- - + + ]
880 : :
881 : 76 : call htetra_q%init(bz2ibz, crystal%gprimd, qlatt, qbz, nqbz, qibz, phdos%nqibz, ierr, errstr, comm)
882 : : !call cwtime_report(" init_tetra", cpu, wall, gflops)
883 [ - + ]: 76 : ABI_CHECK(ierr == 0, errstr)
884 [ - + ]: 76 : ABI_FREE(bz2ibz)
885 : :
886 : : ! Allocate arrays used to store the entire spectrum, Required to calculate tetra weights.
887 : : ! this may change in the future if Matteo refactorizes the tetra weights as sums over k instead of sums over bands
888 [ + - - + : 104400 : ABI_CALLOC(full_phfrq, (3*natom, phdos%nqibz))
+ - - + +
- + - + -
- + - + +
+ + + ]
889 [ - + + - : 369509 : ABI_CALLOC(full_veloc, (3,3*natom, phdos%nqibz))
- + + - +
- + - - +
- + + + +
+ + + ]
890 [ - + + - : 228 : ABI_MALLOC(full_phangmom, (3, 3*natom, phdos%nqibz))
- + + - -
+ - + ]
891 [ - + + - : 380 : ABI_MALLOC_OR_DIE(full_eigvec, (2, 3, natom, 3*natom, phdos%nqibz), ierr)
- + + - -
+ + - + -
+ - + - -
+ - + ]
892 [ + + + + : 1885406 : full_eigvec = zero
+ + + + +
+ ]
893 : : end if ! tetra
894 : :
895 [ - + ]: 177 : ABI_FREE(bz2ibz_smap)
896 [ - + ]: 177 : ABI_FREE(new_shiftq)
897 : :
898 : : ! MPI Sum over irreducible q-points then sync the following integrals:
899 : : ! speedofsound, nsmallq
900 : : ! wminmax and count_wminmax
901 : : ! if gauss: %phdos, %msqd_dos_atom
902 : : ! if tetra: full_phfrq, full_eigvec, full_phangmom, %phdos_int
903 : :
904 : 177 : nsmallq = zero; speedofsound = zero
905 : 177 : wminmax = [huge(one), -huge(one)]; count_wminmax = 0
906 : :
907 [ + + ]: 129342 : do iq_ibz=1,phdos%nqibz
908 [ + - ]: 129165 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
909 : :
910 : : ! Fourier interpolation (keep track of min/max to decide if initial mesh was large enough)
911 : : ! get velocities from this call as well for iq_ibz
912 : : ! the if clause saves a bit of time if we do not need to project anything later
913 [ + - ]: 129165 : if (phdos%n_normal_vec_dmm > 0) then
914 : 129165 : call ifc%fourq(crystal, qibz(:,iq_ibz), phfrq, displ, out_eigvec=eigvec, dwdq=veloc_1q)
915 : : else
916 : 0 : call ifc%fourq(crystal, qibz(:,iq_ibz), phfrq, displ, out_eigvec=eigvec)
917 : : end if
918 : :
919 [ + - - + : 946287 : wminmax(1) = min(wminmax(1), minval(phfrq))
- - + + ]
920 [ + + ]: 129165 : if (wminmax(1) < phdos%omega(1)) count_wminmax(1) = count_wminmax(1) + 1
921 [ + - - + : 946287 : wminmax(2) = max(wminmax(2), maxval(phfrq))
- - + + ]
922 [ + + ]: 129165 : if (wminmax(2) > phdos%omega(nomega)) count_wminmax(2) = count_wminmax(2) + 1
923 : :
924 [ + + ]: 516660 : normq = sum(qibz(:,iq_ibz) ** 2)
925 [ + + + + ]: 129165 : if (normq < max_smallq .and. normq > tol6) then
926 : 10369 : call phdos_calc_vsound(eigvec, crystal%gmet, natom, phfrq, qibz(:,iq_ibz), speedofsound_)
927 [ + + ]: 41476 : speedofsound = speedofsound + speedofsound_ * wtq_ibz(iq_ibz)
928 : 10369 : nsmallq = nsmallq + wtq_ibz(iq_ibz)
929 : : end if
930 : :
931 : 129165 : dos_maxmode_ = 3*natom
932 [ + + ]: 129165 : if (present(dos_maxmode)) then
933 [ - + ]: 1982 : if (dos_maxmode > 0 .and. dos_maxmode < 3*natom) then
934 : 129165 : dos_maxmode_ = dos_maxmode
935 : : end if
936 : : end if
937 : :
938 : : ! Optionally skip the three acoustic modes at Gamma. They are the three modes with the
939 : : ! smallest |frequency|, not necessarily the lowest ones (unstable modes have w<0).
940 [ + + ]: 817122 : skip_mode = .false.
941 [ + + ]: 129165 : if (present(exclude_gamma_acoustic)) then
942 [ + - + + ]: 111458 : if (exclude_gamma_acoustic .and. normq < tol12) then
943 [ + + ]: 752 : abs_phfrq = abs(phfrq)
944 [ + + ]: 368 : do ii=1,3
945 [ + - - + : 2256 : imin = minloc(abs_phfrq, dim=1)
- - + + +
+ ]
946 : 276 : skip_mode(imin) = .true.
947 : 368 : abs_phfrq(imin) = huge(one)
948 : : end do
949 : : end if
950 : : end if
951 : :
952 : 177 : select case (prtdos)
953 : : case (1)
954 [ + + ]: 713102 : do imode=1,dos_maxmode_
955 [ + + ]: 599562 : if (skip_mode(imode)) cycle
956 : : ! Precompute \delta(w - w_{qnu}) * weight(q)
957 [ + - - + : 341796966 : xvals = (phdos%omega(:) - phfrq(imode)) * gaussfactor
- - + + ]
958 [ + + + + ]: 341197680 : where (abs(xvals) < gaussmaxarg)
959 : : gvals_wtq = gaussprefactor * exp(-xvals*xvals) * wtq_ibz(iq_ibz)
960 : : elsewhere
961 : : gvals_wtq = zero
962 : : end where
963 : :
964 : : ! Accumulate PHDOS
965 [ + + ]: 341197680 : phdos%phdos(:) = phdos%phdos(:) + gvals_wtq(:)
966 : :
967 : : ! Rotate e(q) to get e(Sq) to account for symmetrical q-points in BZ.
968 : : ! eigenvectors indeed are not invariant under rotation. See e.g. Eq 39-40 of PhysRevB.76.165108 [[cite:Giustino2007]].
969 : : ! In principle there's a phase due to nonsymmorphic translations but we here need |e(Sq)_iatom|**2
970 [ + + + + ]: 5570094 : syme2_xyza = zero
971 [ + + ]: 1841988 : do iat=1,natom
972 [ + + ]: 18520524 : do isym=1,my_nsym
973 : 16678536 : jat = crystal%indsym(4,isym,iat)
974 : : syme2_xyza(:,jat) = syme2_xyza(:,jat) + &
975 : : matmul(symcart(:,:,isym), eigvec(1,:,iat,imode)) ** 2 + &
976 [ + + + + : 468241710 : matmul(symcart(:,:,isym), eigvec(2,:,iat,imode)) ** 2
+ + + + +
+ ]
977 : : end do
978 : : end do
979 [ + + + + ]: 5570094 : syme2_xyza = syme2_xyza / my_nsym
980 : :
981 : : ! Accumulate PJDOS
982 [ + + ]: 1841988 : do iat=1,natom
983 [ + + ]: 5570094 : do idir=1,3
984 [ + + ]: 2152020501 : phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * gvals_wtq
985 : : end do
986 : : end do
987 : :
988 : : ! Accumulate outer product of displacement vectors
989 : : ! NB: only accumulate real part. e(-q) = e(q)* the full sum over the BZ guarantees Im=0
990 : : ! this sum only does irreducible points: the matrix is symmetrized below
991 : : ! msqd_atom_tmp has units of bohr^2 / Ha as gaussval ~ 1/smear ~ 1/Ha
992 [ + + ]: 1841988 : do iat=1,natom
993 : 1242702 : msqd_atom_tmp = zero
994 [ + + ]: 4970808 : do idir=1,3
995 [ + + ]: 16155126 : do jdir=1,3
996 : : msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( &
997 : : eigvec(1,idir,iat,imode) * eigvec(1,jdir,iat,imode) &
998 : 14912424 : + eigvec(2,idir,iat,imode) * eigvec(2,jdir,iat,imode) )
999 : : end do
1000 : : end do
1001 : :
1002 : : ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ.
1003 : : ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately
1004 : : ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause
1005 : : ! NB: looks consistent with the sym in harmonic thermo, just used in opposite
1006 : : ! direction for symops: symrel here instead of symrec and the inverse of indsym in harmonic_thermo
1007 [ + + ]: 18520524 : do isym=1, my_nsym
1008 : : !temp_33 = matmul( (symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))) )
1009 : : ! MG Version
1010 [ + + + + : 1517746776 : temp_33 = matmul( (transpose(symcart(:,:,isym))), matmul(msqd_atom_tmp, symcart(:,:,isym)) )
+ + + + +
+ + + + +
+ + ]
1011 [ + + + + ]: 216820968 : temp_33 = temp_33 / my_nsym
1012 : 16678536 : jat = crystal%indsym(4,isym,iat)
1013 [ + + ]: 67956846 : do idir=1,3
1014 [ + + ]: 216820968 : do jdir=1,3
1015 : : phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + &
1016 [ + + ]: >10561*10^7 : temp_33(idir, jdir) * gvals_wtq
1017 : : end do
1018 : : end do
1019 : : end do
1020 : :
1021 : : end do ! iat msqd loop
1022 : :
1023 : : ! DOS including velocity projected onto interface normal vector
1024 [ + + ]: 2397144 : veloc = veloc_1q(:,imode)
1025 [ + + ]: 4907828 : do icart = 1, phdos%n_normal_vec_dmm
1026 : : ! projection of the velocity onto chosen normal vector for the interface in DMM.
1027 : : ! NB: The sign can change, but norm of scalar product is the correct recipe.
1028 : : projfact = zero
1029 : : ! as we are summing only over irred q-points, need to add all contributions
1030 : : ! from the velocities on the star of q_ibz
1031 [ + + ]: 71059422 : do isym=1,my_nsym
1032 [ + + + + : 1074025722 : projfact = projfact + abs(dot_product(matmul(symcart(:,:,isym),veloc),phdos%normal_vec_dmm(:,icart)))
+ + ]
1033 : : end do
1034 [ + + ]: 2388983322 : phdos%phdos_dmm(:,icart) = phdos%phdos_dmm(:,icart) + projfact * gvals_wtq(:) / my_nsym
1035 : : end do
1036 : :
1037 : : end do ! imode
1038 : :
1039 : : case (2)
1040 : : ! Tetrahedra; Save phonon frequencies, eigenvectors and angular momentum.
1041 : : ! Sum is done after the loops over the two meshes.
1042 [ + + ]: 104020 : full_phfrq(:,iq_ibz) = phfrq(:)
1043 [ + + + + ]: 369205 : full_veloc(:,:,iq_ibz) = veloc_1q(:,:)
1044 [ + + + + : 1885330 : full_eigvec(:,:,:,:,iq_ibz) = eigvec
+ + + + ]
1045 [ + + + + ]: 369205 : full_phangmom(:,:,iq_ibz) = phangmom
1046 : :
1047 : : case default
1048 [ + + - ]: 129165 : ABI_ERROR(sjoin("Wrong value for prtdos:", itoa(prtdos)))
1049 : : end select
1050 : : end do ! iq_ibz
1051 : :
1052 [ - + ]: 177 : ABI_FREE(qbz)
1053 [ - + ]: 177 : ABI_FREE(gvals_wtq)
1054 [ - + ]: 177 : ABI_FREE(xvals)
1055 : :
1056 : 177 : call xmpi_sum_master(nsmallq, master, comm, ierr)
1057 : 177 : call xmpi_sum_master(speedofsound, master, comm, ierr)
1058 : :
1059 : : !call cwtime_report(" phdos", cpu, wall, gflops)
1060 : :
1061 [ + - ]: 177 : if (my_rank == master) then
1062 [ + + ]: 177 : if (nsmallq > tol10) then
1063 : : ! Write info about speed of sound
1064 [ + + ]: 700 : speedofsound = speedofsound / nsmallq
1065 : : write (msg,'(a,E20.10,3a,F16.4,2a)') &
1066 [ + + ]: 700 : ' Average speed of sound partial sums: ', third*sum(speedofsound), ' (at units)',ch10, &
1067 [ + + ]: 875 : '- = ', third*sum(speedofsound) * Bohr_Ang * 1.d-13 / Time_Sec, ' [km/s]',ch10
1068 [ + + ]: 175 : if (prtout_) call wrtout(units, msg)
1069 : :
1070 : : ! Debye frequency = vs * (6 pi^2 natom / ucvol)**1/3
1071 [ + + ]: 700 : debyefreq = third*sum(speedofsound) * (six*pi**2/crystal%ucvol)**(1./3.)
1072 : : write (msg,'(a,E20.10,3a,E20.10,a)') &
1073 : 175 : ' Debye frequency from partial sums: ', debyefreq, ' (Ha)',ch10, &
1074 : 350 : '- = ', debyefreq*Ha_THz, ' (THz)'
1075 [ + + ]: 175 : if (prtout_) call wrtout(units, msg)
1076 : :
1077 : : ! Debye temperature = hbar * Debye frequency / kb
1078 : 175 : write (msg,'(a,E20.10,2a)') '-Debye temperature from partial sums: ', debyefreq*Ha_K, ' (K)', ch10
1079 [ + + ]: 175 : if (prtout_) call wrtout(units, msg)
1080 : : end if
1081 : : end if
1082 : :
1083 [ + + ]: 177 : if (prtdos == 2) then
1084 : 76 : call cwtime(cpu, wall, gflops, "start")
1085 : : ! Finalize integration with tetrahedra
1086 : : ! All the data are contained in full_phfrq, full_eigvec and full_phangmom.
1087 : 76 : call xmpi_sum(full_phfrq, comm, ierr)
1088 : 76 : call xmpi_sum(full_veloc, comm, ierr)
1089 : 76 : call xmpi_sum(full_eigvec, comm, ierr)
1090 : 76 : call xmpi_sum(full_phangmom, comm, ierr)
1091 : :
1092 [ - + + - : 228 : ABI_MALLOC(tmp_phfrq, (phdos%nqibz))
+ - + - -
+ - + ]
1093 : :
1094 : 0 : do
1095 [ - + + - : 228 : ABI_MALLOC(wdt, (phdos%nomega, 2))
+ - + - -
+ - + ]
1096 [ - + + - : 228 : ABI_MALLOC(energies, (phdos%nomega))
+ - + - -
+ - + -
+ ]
1097 [ - + ]: 76 : energies = linspace(phdos%omega_min, phdos%omega_max, phdos%nomega)
1098 : :
1099 [ + + ]: 15701 : do iq_ibz=1,phdos%nqibz
1100 [ + - ]: 15625 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism
1101 : :
1102 [ + + ]: 104096 : do imode=1,dos_maxmode_
1103 : : ! Compute the weights for this q-point using tetrahedron
1104 [ + + ]: 20158392 : tmp_phfrq(:) = full_phfrq(imode,:)
1105 : 88395 : call htetra_q%get_onewk_wvals(iq_ibz,bcorr0,phdos%nomega,energies,max_occ1,phdos%nqibz,tmp_phfrq,wdt)
1106 [ + + + + ]: 161295981 : wdt = wdt * wtq_ibz(iq_ibz)
1107 : :
1108 : : ! Accumulate DOS/IDOS
1109 [ + + ]: 80603793 : phdos%phdos(:) = phdos%phdos(:) + wdt(:, 1)
1110 [ + + ]: 80603793 : phdos%phdos_int(:) = phdos%phdos_int(:) + wdt(:, 2)
1111 : :
1112 : : ! Rotate e(q) to get e(Sq) to account for other q-points in BZ. See notes in gaussian branch
1113 [ + + + + ]: 800919 : syme2_xyza = zero
1114 [ + + ]: 266526 : do iat=1,natom
1115 [ + + ]: 7665588 : do isym=1,my_nsym
1116 : 7399062 : jat = crystal%indsym(4,isym,iat)
1117 : : syme2_xyza(:,jat) = syme2_xyza(:,jat) + &
1118 : : matmul(symcart(:,:,isym), full_eigvec(1,:,iat,imode,iq_ibz)) ** 2 + &
1119 [ + + + + : 207351867 : matmul(symcart(:,:,isym), full_eigvec(2,:,iat,imode,iq_ibz)) ** 2
+ - + + +
+ + - +
+ ]
1120 : : end do
1121 : : end do
1122 [ + + + + ]: 800919 : syme2_xyza = syme2_xyza / my_nsym
1123 : :
1124 [ + + ]: 266526 : do iat=1,natom
1125 [ + + ]: 800919 : do idir=1,3
1126 [ + + ]: 597389409 : phdos%pjdos(:,idir,iat) = phdos%pjdos(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,1)
1127 [ + + ]: 597567540 : phdos%pjdos_int(:,idir,iat) = phdos%pjdos_int(:,idir,iat) + syme2_xyza(idir,iat) * wdt(:,2)
1128 : : end do
1129 : : end do
1130 : :
1131 [ + + ]: 266526 : do iat=1,natom
1132 : :
1133 : : ! Accumulate outer product of displacement vectors
1134 : 178131 : msqd_atom_tmp = zero
1135 [ + + ]: 712524 : do idir=1,3
1136 [ + + ]: 2315703 : do jdir=1,3
1137 : : msqd_atom_tmp(jdir,idir) = msqd_atom_tmp(jdir,idir) + ( &
1138 : : full_eigvec(1,idir,iat,imode,iq_ibz) * full_eigvec(1,jdir,iat,imode,iq_ibz) &
1139 : 2137572 : + full_eigvec(2,idir,iat,imode,iq_ibz) * full_eigvec(2,jdir,iat,imode,iq_ibz) )
1140 : : end do
1141 : : end do
1142 : :
1143 : : ! Symmetrize matrices to get full sum of tensor over all BZ, not just IBZ.
1144 : : ! the atom is not necessarily invariant under symops, so these contributions should be added to each iat separately
1145 : : ! normalization by nsym is done at the end outside the iqpt loop and after the tetrahedron clause
1146 : : ! from loops above only the eigvec are kept and not the displ, so we still have to divide by the masses
1147 : : ! TODO: need to check the direction of the symcart vs transpose or inverse, given that jat is the pre-image of iat...
1148 [ + + ]: 7665588 : do isym=1,my_nsym
1149 : : !temp_33 = matmul((symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))))
1150 : : ! MG Version
1151 [ + + + + : 673314642 : temp_33 = matmul( (transpose(symcart(:,:,isym))), matmul(msqd_atom_tmp, symcart(:,:,isym)) )
+ + + + +
+ + + + +
+ + ]
1152 [ + + + + ]: 96187806 : temp_33 = temp_33 / my_nsym
1153 : 7399062 : jat = crystal%indsym(4,isym,iat)
1154 [ + + ]: 29774379 : do idir=1,3
1155 [ + + ]: 96187806 : do jdir=1,3
1156 : : phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + &
1157 [ + + ]: 61356659490 : temp_33(idir, jdir) * wdt(:,1)
1158 : : end do
1159 : : end do
1160 : : end do
1161 : : end do ! iat msqd accumulation
1162 : :
1163 : : ! DOS including velocity projected onto interface normal vector
1164 [ + + ]: 353580 : veloc = full_veloc(:,imode,iq_ibz)
1165 [ + + ]: 722785 : do icart = 1, phdos%n_normal_vec_dmm
1166 : : ! projection of the velocity onto chosen normal vector for the interface in DMM.
1167 : : ! NB: The sign can change, but norm of scalar product is the correct recipe.
1168 : : projfact = zero
1169 : : ! as we are summing only over irred q-points, need to add all contributions
1170 : : ! from the velocities on the star of q_ibz
1171 [ + + ]: 27404727 : do isym=1,my_nsym
1172 [ + + + + : 429194157 : projfact = projfact + abs(dot_product(matmul(symcart(:,:,isym),veloc),phdos%normal_vec_dmm(:,icart)))
+ + ]
1173 : : end do
1174 [ + + ]: 564314946 : phdos%phdos_dmm(:,icart) = phdos%phdos_dmm(:,icart) + projfact * wdt(:,1) / my_nsym
1175 : : end do
1176 : :
1177 : : end do ! imode
1178 : : end do ! iq_ibz
1179 : :
1180 [ - + ]: 76 : if (refine_dosdeltae) then
1181 : : ! HM: Check if the integration of the DOS is correct, otherwise half dos%deltae and re-run
1182 : : ! MG FIXME: This won't work in parallel because we still have to call xmpi_sum
1183 : 0 : call ctrap(phdos%nomega, phdos%phdos, phdos%omega_step, phdos_int)
1184 [ # # ]: 0 : if (abs(phdos_int - crystal%natom*3) > tol2) then
1185 : 0 : write(msg,'(a,f6.2,a,i4,2a,e10.3,a,e10.3)') "The value of the integral is", phdos_int, &
1186 : 0 : " but it should be", crystal%natom*3, ch10,&
1187 : 0 : "I will decrease dosdeltae from: ", dosdeltae, " to: ", dosdeltae/two
1188 : 0 : ABI_WARNING(msg)
1189 : 0 : dosdeltae = dosdeltae / two
1190 : 0 : call phdos%free()
1191 : 0 : call phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos)
1192 : 0 : nomega = phdos%nomega
1193 [ # # ]: 0 : ABI_FREE(wdt)
1194 [ # # ]: 0 : ABI_FREE(energies)
1195 : : cycle
1196 : : endif
1197 : : endif
1198 : :
1199 : : exit
1200 : : end do
1201 : :
1202 [ - + ]: 76 : ABI_FREE(energies)
1203 [ - + ]: 76 : ABI_FREE(wdt)
1204 : :
1205 : : ! Make eigvec into phonon displacements.
1206 [ + + ]: 225 : do iat = 1, natom
1207 [ + + + + : 1811000 : full_eigvec(:,:,iat,:,:) = full_eigvec(:,:,iat,:,:) / sqrt(phdos%atom_mass(iat))
+ + + + ]
1208 : : end do
1209 : :
1210 : 76 : call cwtime_report(" tetra accumulate", cpu, wall, gflops)
1211 : :
1212 [ + - - + ]: 76 : if (my_rank == master .and. len_trim(prefix) > 0) then
1213 [ # # ]: 0 : NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHIBZ.nc"), xmpi_comm_self), "Creating PHIBZ")
1214 [ # # ]: 0 : NCF_CHECK(crystal%ncwrite(ncid))
1215 : 0 : call phonons_ncwrite(ncid, natom, phdos%nqibz, qibz, wtq_ibz, full_phfrq, full_eigvec, full_phangmom)
1216 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
1217 : : end if
1218 : :
1219 : : ! Immediately free this - it contains displ and not eigvec at this stage
1220 [ - + ]: 76 : ABI_FREE(full_eigvec)
1221 [ - + ]: 76 : ABI_FREE(full_phfrq)
1222 [ - + ]: 76 : ABI_FREE(full_veloc)
1223 [ - + ]: 76 : ABI_FREE(full_phangmom)
1224 [ - + ]: 76 : ABI_FREE(tmp_phfrq)
1225 : 228 : call htetra_q%free()
1226 : : else
1227 : 101 : ABI_WARNING('The netcdf PHIBZ file is only output for tetrahedron integration and DOS calculations')
1228 : : end if ! tetrahedra
1229 : :
1230 : : ! Test if the initial mesh was large enough
1231 : 177 : call xmpi_sum(count_wminmax, comm, ierr)
1232 : 354 : call xmpi_min(wminmax(1), rtmp, comm, ierr); wminmax(1) = rtmp
1233 : 354 : call xmpi_max(wminmax(2), rtmp, comm, ierr); wminmax(2) = rtmp
1234 : :
1235 : 177 : call xmpi_sum(phdos%phdos, comm, ierr)
1236 : 177 : call xmpi_sum(phdos%phdos_dmm, comm, ierr)
1237 : 177 : call xmpi_sum(phdos%msqd_dos_atom, comm, ierr)
1238 : 177 : call xmpi_sum(phdos%pjdos, comm, ierr)
1239 [ + + ]: 177 : if (prtdos == 2) then
1240 : : ! Reduce integrals if tetra, gauss will compute idos with simpson.
1241 : 76 : call xmpi_sum(phdos%phdos_int, comm, ierr)
1242 : 76 : call xmpi_sum(phdos%pjdos_int, comm, ierr)
1243 : : end if
1244 : :
1245 : : ! for dmm phdos need to normalize by q point integration element
1246 [ + + + + ]: 1438964 : phdos%phdos_dmm = phdos%phdos_dmm * two_pi**3 / crystal%ucvol
1247 : :
1248 : : ! ===============================
1249 : : ! === Compute Integrated PDOS ===
1250 : : ! ===============================
1251 [ - + + - : 903813 : ABI_CALLOC(phdos%pjdos_rc_type, (nomega, 3, crystal%ntypat))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
1252 [ - + + - : 301661 : ABI_CALLOC(phdos%pjdos_type, (nomega, crystal%ntypat))
- + + - +
- + - - +
- + - + +
+ + + ]
1253 [ - + + - : 301484 : ABI_CALLOC(phdos%pjdos_type_int, (nomega, crystal%ntypat))
- + + - -
+ - + - +
+ + + + ]
1254 : :
1255 [ + + ]: 555 : do iat=1,natom
1256 : 378 : itype = crystal%typat(iat)
1257 [ + + ]: 480477 : do io=1,phdos%nomega
1258 [ + + ]: 1920396 : phdos%pjdos_rc_type(io,:,itype) = phdos%pjdos_rc_type(io,:,itype) + phdos%pjdos(io,:,iat)
1259 [ + + ]: 1920774 : phdos%pjdos_type(io,itype) = phdos%pjdos_type(io,itype) + sum(phdos%pjdos(io,:,iat))
1260 : : end do
1261 [ + + ]: 555 : if (prtdos == 2) then
1262 [ + + ]: 199348 : do io=1,phdos%nomega
1263 [ + + ]: 796945 : phdos%pjdos_type_int(io,itype) = phdos%pjdos_type_int(io,itype) + sum(phdos%pjdos_int(io,:,iat))
1264 : : end do
1265 : : end if
1266 : : end do
1267 : :
1268 : : ! Evaluate IDOS using simple simpson integration
1269 : : ! In principle one could use derf.F90, just to be consistent ...
1270 [ + + ]: 177 : if (prtdos == 1) then
1271 : 101 : call simpson_int(phdos%nomega,phdos%omega_step, phdos%phdos, phdos%phdos_int)
1272 [ + + ]: 330 : do iat=1,natom
1273 [ + + ]: 1017 : do idir=1,3
1274 : 916 : call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos(:,idir,iat), phdos%pjdos_int(:,idir,iat))
1275 : : end do
1276 : : end do
1277 [ + + ]: 224 : do itype=1,crystal%ntypat
1278 : 224 : call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos_type(:,itype), phdos%pjdos_type_int(:,itype))
1279 : : end do
1280 : : end if
1281 : :
1282 [ - + ]: 177 : ABI_FREE(qibz)
1283 [ - + ]: 177 : ABI_FREE(wtq_ibz)
1284 : :
1285 : 177 : call cwtime_report(" phdos_init", cpu_all, wall_all, gflops_all)
1286 : : DBG_EXIT("COLL")
1287 : :
1288 [ - + - + : 531 : end subroutine phdos_init
- + - + -
- - - - +
- - - - -
+ - + - +
- + - + ]
1289 : : !!***
1290 : :
1291 : : !----------------------------------------------------------------------
1292 : :
1293 : : !!****f* m_phonons/zacharias_supercell_make
1294 : : !! NAME
1295 : : !! zacharias_supercell_make
1296 : : !!
1297 : : !! FUNCTION
1298 : : !! Construct an optimally thermalized supercell following Zacharias and Giustino
1299 : : !! See: PRB 94 075125 (2016) [[cite:Zacharias2016]]
1300 : : !!
1301 : : !! INPUTS
1302 : : !!
1303 : : !! OUTPUT
1304 : : !!
1305 : : !! SOURCE
1306 : :
1307 [ + + - + : 21 : subroutine zacharias_supercell_make(Crystal, Ifc, ntemper, rlatt, tempermin, temperinc, thm_scells)
- + - + -
+ - + - +
- + + + +
- ]
1308 : :
1309 : : !Arguments ------------------------------------
1310 : : !scalars
1311 : : integer, intent(in) :: ntemper
1312 : : integer, intent(in) :: rlatt(3,3)
1313 : : real(dp), intent(in) :: tempermin, temperinc
1314 : : type(crystal_t),intent(in) :: Crystal
1315 : : type(ifc_type),intent(in) :: Ifc
1316 : : type(supercell_type), intent(out) :: thm_scells(ntemper)
1317 : :
1318 : : !Local variables-------------------------------
1319 : : !scalars
1320 : : integer :: iq, nqibz, nqbz, qptopt1, imode, itemper, ierr, jmode
1321 : : real(dp) :: temperature_K, temperature, modesign, sigma, freeze_displ
1322 : : !arrays
1323 : 1 : integer, allocatable :: modeindex(:)
1324 : 1 : real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt
1325 : 1 : real(dp), allocatable :: qbz(:,:), qibz(:,:), wtq_ibz(:)
1326 : 1 : real(dp), allocatable :: phfrq_allq(:), phdispl_allq(:,:,:,:,:)
1327 : 1 : real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:)
1328 : 1 : real(dp), allocatable :: phdispl1(:,:,:)
1329 : : character (len=500) :: msg
1330 : : ! *************************************************************************
1331 : :
1332 : : ! check inputs
1333 : : ! TODO: add check that all rlatt are the same on input
1334 : :
1335 : : if (rlatt(1,2)/=0 .or. rlatt(1,3)/=0 .or. rlatt(2,3)/=0 .or. &
1336 [ + - + - : 1 : rlatt(2,1)/=0 .or. rlatt(3,1)/=0 .or. rlatt(3,2)/=0) then
+ - + - +
- - + ]
1337 : 0 : write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', &
1338 : 0 : ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning '
1339 : 0 : ABI_WARNING(msg)
1340 : 0 : return
1341 : : end if
1342 : :
1343 : : ! build qpoint grid used for the Fourier interpolation (use no syms for the moment!)
1344 : 1 : qptopt1 = 3
1345 : :
1346 : : ! for the moment do not allow shifted q grids. We are interpolating anyway, so it will always work
1347 [ - + ]: 1 : ABI_MALLOC(qshft,(3,1))
1348 [ + + ]: 4 : qshft(:,1)=zero
1349 : :
1350 : : ! This call will set nqibz, IBZ and BZ arrays
1351 : : call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, &
1352 : 1 : nqibz, qibz, wtq_ibz, nqbz, qbz) ! new_kptrlatt, new_shiftk) ! Optional
1353 [ - + ]: 1 : ABI_FREE(qshft)
1354 : :
1355 : : ! allocate arrays with all of the q, omega, and displacement vectors
1356 [ + - - + : 3 : ABI_MALLOC_OR_DIE(phfrq_allq, (3*Crystal%natom*nqibz), ierr)
+ - + - +
- + - - +
- + ]
1357 [ - + + - : 5 : ABI_MALLOC_OR_DIE(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr)
+ - - + +
- - + + -
+ - + - +
- - + -
+ ]
1358 : :
1359 [ - + + - : 3 : ABI_MALLOC(phfrq, (3*Crystal%natom))
+ - + - -
+ - + ]
1360 [ - + + - : 4 : ABI_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom))
- + + - +
- + - - +
- + ]
1361 [ - + + - : 3 : ABI_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom))
- + + - -
+ - + ]
1362 : :
1363 : : ! loop over q to get all frequencies and displacement vectors
1364 [ - + + - : 3 : ABI_MALLOC(modeindex, (nqibz*3*Crystal%natom))
+ - + - -
+ - + ]
1365 : 1 : imode = 0
1366 [ + + ]: 9 : do iq = 1, nqibz
1367 : : ! Fourier interpolation.
1368 : 8 : call ifc%fourq(Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec)
1369 [ + + ]: 152 : phfrq_allq((iq-1)*3*Crystal%natom+1 : iq*3*Crystal%natom) = phfrq
1370 [ + + + + : 8792 : phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl
+ + + + ]
1371 [ + + ]: 153 : do jmode = 1, 3*Crystal%natom
1372 : 144 : imode = imode + 1
1373 : 152 : modeindex(imode) = imode
1374 : : end do
1375 : : end do
1376 [ - + ]: 1 : ABI_FREE(phfrq)
1377 [ - + ]: 1 : ABI_FREE(pheigvec)
1378 [ - + ]: 1 : ABI_FREE(phdispl)
1379 : :
1380 : : ! sort modes in whole list: get indirect indexing for qbz and displ
1381 : 1 : call sort_dp(nqibz*3*Crystal%natom, phfrq_allq, modeindex, tol10)
1382 : : ! NB: phfrq is sorted now, but displ and qibz will have to be indexed indirectly with modeindex
1383 : :
1384 : : ! only diagonal supercell case for the moment
1385 [ + + ]: 11 : do itemper = 1, ntemper
1386 : 11 : call thm_scells(itemper)%init(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl)
1387 : : end do
1388 : :
1389 : : ! precalculate phase factors???
1390 [ - + + - : 3 : ABI_MALLOC(phdispl1, (2, 3, Crystal%natom))
+ - + - -
+ - + ]
1391 : : ! for all modes at all q in whole list, sorted
1392 : 145 : modesign=one
1393 : :
1394 [ + + ]: 145 : do imode = 1, 3*Crystal%natom*nqibz
1395 : : ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma)
1396 [ + + ]: 144 : if (phfrq_allq(imode) < tol10) cycle
1397 : :
1398 [ + + ]: 141 : iq = ceiling(dble(modeindex(imode))/dble(3*Crystal%natom))
1399 : 141 : jmode = modeindex(imode) - (iq-1)*3*Crystal%natom
1400 [ - + + + : 8742 : phdispl1 = phdispl_allq(:,:,:,jmode,iq)
+ + + + ]
1401 : :
1402 : : ! loop over temperatures
1403 [ + + ]: 1551 : do itemper = 1, ntemper
1404 : 1410 : temperature_K = tempermin + dble(itemper-1)*temperinc ! this is in Kelvin
1405 : 1410 : temperature = temperature_K / Ha_K
1406 : :
1407 : : ! trick supercell object into using present q point
1408 [ + + ]: 5640 : thm_scells(itemper)%qphon(:) = qibz(:,iq)
1409 : :
1410 : : ! find thermal displacement amplitude eq 4 of Zacharias
1411 : : ! combined with l_nu,q expression in paragraph before
1412 : 1410 : sigma = sqrt((bose_einstein(phfrq_allq(imode), temperature) + half)/phfrq_allq(imode))
1413 : :
1414 : : ! add displacement for this mode to supercell positions eq 5 of Zacharias
1415 : 1410 : freeze_displ = modesign * sigma
1416 : 1551 : call thm_scells(itemper)%freeze_displ(phdispl1(:,:,:), freeze_displ)
1417 : : end do !itemper
1418 : :
1419 : : ! this is the prescription: flip sign for each successive mode in full
1420 : : ! spectrum, to cancel electron phonon coupling to 1st order
1421 : : ! (hopefully 3rd order as well)
1422 : 145 : modesign = -modesign
1423 : : end do !imode
1424 : :
1425 [ - + ]: 1 : ABI_FREE(modeindex)
1426 [ - + ]: 1 : ABI_FREE(phfrq_allq)
1427 [ - + ]: 1 : ABI_FREE(phdispl_allq)
1428 [ - + ]: 1 : ABI_FREE(phdispl1)
1429 [ - + ]: 1 : ABI_FREE(qibz)
1430 [ - + ]: 1 : ABI_FREE(qbz)
1431 [ - + ]: 1 : ABI_FREE(wtq_ibz)
1432 : :
1433 [ - + - + : 1 : end subroutine zacharias_supercell_make
- + ]
1434 : : !!***
1435 : :
1436 : : !----------------------------------------------------------------------
1437 : :
1438 : : !!****f* m_phonons/thermal_supercell_make
1439 : : !! NAME
1440 : : !! thermal_supercell_make
1441 : : !!
1442 : : !! FUNCTION
1443 : : !! Construct a random thermalized supercell configuration, as in TDEP
1444 : : !! main function is for training set generation in multibinit
1445 : : !!
1446 : : !! INPUTS
1447 : : !! Crystal = crystal object with rprim etc...
1448 : : !! Ifc = interatomic force constants object from anaddb
1449 : : !! option = option to deal with negative frequency -> Bose factor explodes (eg acoustic at Gamma)
1450 : : !! several philosophies to be implemented for the unstable modes:
1451 : : !! option == 1 => ignore
1452 : : !! option == 2 => populate them according to a default amplitude
1453 : : !! option == 3 => populate according to their modulus squared
1454 : : !! option == 4 => USER defined value(s), require namplitude and amplitude
1455 : : !! nconfig = number of requested configurations
1456 : : !! rlatt = matrix of conversion for supercell (3 0 0 0 3 0 0 0 3 for example)
1457 : : !! temperature_K = temperature in Kelvin
1458 : : !! nqpt = number of q-point
1459 : : !! namplitude = number of amplitude provided by the user
1460 : : !! amplitudes(namplitude) = list of the amplitudes of the unstable phonons
1461 : : !! amplitudes(1:3,iamplitude) = qpt
1462 : : !! amplitudes(4,iamplitude) = mode
1463 : : !! amplitudes(5,iamplitude) = amplitude
1464 : : !!
1465 : : !! OUTPUT
1466 : : !! thm_scells = array of configurations with thermalized supercells
1467 : : !!
1468 : : !! SOURCE
1469 : :
1470 [ # # # # : 0 : subroutine thermal_supercell_make(amplitudes,Crystal, Ifc,namplitude, nconfig,option,&
# # # # #
# # # # #
# # # # #
# ]
1471 : 0 : rlatt, temperature_K, thm_scells)
1472 : :
1473 : : !Arguments ------------------------------------
1474 : : !scalars
1475 : : integer, intent(in) :: option,nconfig
1476 : : integer, intent(in) :: rlatt(3,3)
1477 : : real(dp), intent(in) :: temperature_K
1478 : : type(crystal_t),intent(in) :: Crystal
1479 : : type(ifc_type),intent(in) :: Ifc
1480 : : type(supercell_type), intent(out) :: thm_scells(nconfig)
1481 : : integer,intent(in) :: namplitude
1482 : : !Local variables-------------------------------
1483 : : !scalars
1484 : : integer :: iq, nqibz, nqbz, qptopt1, iampl ,imode, ierr, iconfig
1485 : : real(dp) :: temperature, sigma, freeze_displ
1486 : : real(dp) :: rand !, rand1, rand2
1487 : : real(dp),intent(in):: amplitudes(5,namplitude)
1488 : : !arrays
1489 : 0 : real(dp), allocatable :: qshft(:,:) ! dummy with 2 dimensions for call to kpts_ibz_from_kptrlatt
1490 : 0 : real(dp), allocatable :: qbz(:,:), qibz(:,:), wtqibz(:)
1491 : 0 : real(dp), allocatable :: phfrq_allq(:,:), phdispl_allq(:,:,:,:,:)
1492 : 0 : real(dp), allocatable :: phfrq(:), phdispl(:,:,:,:),pheigvec(:,:,:,:)
1493 : 0 : real(dp), allocatable :: phdispl1(:,:,:)
1494 : : character (len=500) :: msg
1495 : : ! *************************************************************************
1496 : :
1497 : : ! check inputs
1498 : : ! TODO: add check that all rlatt are the same on input
1499 : : if (rlatt(1,2)/=0 .or. rlatt(1,3)/=0 .or. rlatt(2,3)/=0 .or. &
1500 [ # # # # : 0 : rlatt(2,1)/=0 .or. rlatt(3,1)/=0 .or. rlatt(3,2)/=0) then
# # # # #
# # # ]
1501 : 0 : write (msg, '(4a, 9I6, a)') ' for the moment I have not implemented ', &
1502 : 0 : ' non diagonal supercells.',ch10,' rlatt for temp 1 = ', rlatt, ' Returning '
1503 : 0 : ABI_WARNING(msg)
1504 : 0 : return
1505 : : end if
1506 : :
1507 : 0 : temperature = temperature_K / Ha_K
1508 : :
1509 : : ! build qpoint grid used for the Fourier interpolation.
1510 : : !(use no syms for the moment!)
1511 : 0 : qptopt1 = 3
1512 : :
1513 : : ! for the moment do not allow shifted q grids.
1514 : : ! We are interpolating anyway, so it will always work
1515 [ # # ]: 0 : ABI_MALLOC(qshft,(3,1))
1516 [ # # ]: 0 : qshft(:,1)=zero
1517 : :
1518 : : ! This call will set nqibz, IBZ and BZ arrays
1519 : : call kpts_ibz_from_kptrlatt(crystal, rlatt, qptopt1, 1, qshft, &
1520 : 0 : nqibz, qibz, wtqibz, nqbz, qbz) ! new_kptrlatt, new_shiftk) ! Optional
1521 [ # # ]: 0 : ABI_FREE(qshft)
1522 : :
1523 : : ! allocate arrays wzith all of the q, omega, and displacement vectors
1524 [ # # # # : 0 : ABI_MALLOC_OR_DIE(phfrq_allq, (3*Crystal%natom, nqibz), ierr)
# # # # #
# # # # #
# # # # #
# ]
1525 [ # # # # : 0 : ABI_MALLOC_OR_DIE(phdispl_allq, (2, 3, Crystal%natom, 3*Crystal%natom, nqibz), ierr)
# # # # #
# # # # #
# # # # #
# # # ]
1526 : :
1527 [ # # # # : 0 : ABI_MALLOC(phfrq, (3*Crystal%natom))
# # # # #
# # # ]
1528 [ # # # # : 0 : ABI_MALLOC(phdispl, (2, 3, Crystal%natom, 3*Crystal%natom))
# # # # #
# # # # #
# # ]
1529 [ # # # # : 0 : ABI_MALLOC(pheigvec, (2, 3, Crystal%natom, 3*Crystal%natom))
# # # # #
# # # ]
1530 : :
1531 : : ! loop over q to get all frequencies and displacement vectors
1532 : 0 : imode = 0
1533 [ # # ]: 0 : do iq = 1, nqibz
1534 : : ! Fourier interpolation.
1535 : 0 : call ifc%fourq(Crystal, qibz(:,iq), phfrq, phdispl, out_eigvec=pheigvec)
1536 [ # # ]: 0 : phfrq_allq(1:3*Crystal%natom, iq) = phfrq
1537 [ # # # # : 0 : phdispl_allq(1:2, 1:3, 1:Crystal%natom, 1:3*Crystal%natom, iq) = phdispl
# # # # ]
1538 : : end do
1539 [ # # ]: 0 : ABI_FREE(phfrq)
1540 [ # # ]: 0 : ABI_FREE(pheigvec)
1541 [ # # ]: 0 : ABI_FREE(phdispl)
1542 : :
1543 : : ! only diagonal supercell case for the moment
1544 [ # # ]: 0 : do iconfig = 1, nconfig
1545 : 0 : call thm_scells(iconfig)%init(Crystal%natom, rlatt, Crystal%rprimd, Crystal%typat, Crystal%xcart, Crystal%znucl)
1546 : : end do
1547 : :
1548 : : ! precalculate phase factors???
1549 : :
1550 [ # # # # : 0 : ABI_MALLOC(phdispl1, (2, 3, Crystal%natom))
# # # # #
# # # ]
1551 : :
1552 : : ! for all modes at all q in whole list, sorted
1553 [ # # ]: 0 : do iq = 1, nqibz
1554 [ # # ]: 0 : do imode = 1, 3*Crystal%natom
1555 : :
1556 : : ! skip modes with too low or negative frequency -> Bose factor explodes (eg acoustic at Gamma)
1557 : : ! TODO: check the convergence wrt the tolerance
1558 : : ! several philosophies to be implemented for the unstable modes:
1559 : : ! 1) ignore
1560 : : ! 2) populate them according to a default amplitude
1561 : : ! 3) populate according to their modulus squared
1562 [ # # ]: 0 : if (abs(phfrq_allq(imode, iq))<tol6) cycle
1563 : :
1564 [ # # # # : 0 : phdispl1 = phdispl_allq(:,:,:,imode,iq)
# # # # ]
1565 : :
1566 : : ! loop over configurations
1567 [ # # ]: 0 : do iconfig = 1, nconfig
1568 : :
1569 : : ! trick supercell object into using present q point
1570 [ # # ]: 0 : thm_scells(iconfig)%qphon(:) = qibz(:,iq)
1571 : :
1572 : : ! find thermal displacement amplitude eq 4 of Zacharias
1573 : : ! combined with l_nu,q expression in paragraph before
1574 [ # # ]: 0 : if (phfrq_allq(imode, iq) > tol6) then
1575 : 0 : sigma = sqrt( (bose_einstein(phfrq_allq(imode,iq), temperature) + half)/phfrq_allq(imode,iq))
1576 : : else
1577 : : !Treat negative frequencies
1578 : 0 : select case (option)
1579 : : case (1)
1580 : : !Do not populate
1581 : 0 : sigma = 0._dp
1582 : : case (2)
1583 : : !Default amplitude for all the frequencies
1584 : 0 : sigma = 100._dp
1585 : : case (3)
1586 : : !Absolute value of the frequencies
1587 : 0 : sigma=sqrt((bose_einstein(abs(phfrq_allq(imode,iq)),temperature)+half) / abs(phfrq_allq(imode,iq)))
1588 : : case (4)
1589 : : sigma = 0._dp
1590 : : !Search if the amplitude of this unstable phonon is in the input argument amplitudes
1591 [ # # ]: 0 : do iampl=1,namplitude
1592 : : if(abs(thm_scells(iconfig)%qphon(1) - amplitudes(1,iampl)) < tol8.and.&
1593 : : abs(thm_scells(iconfig)%qphon(2) - amplitudes(2,iampl)) < tol8.and.&
1594 [ # # # # : 0 : abs(thm_scells(iconfig)%qphon(3) - amplitudes(3,iampl)) < tol8.and.&
# # # # ]
1595 : 0 : abs(imode - amplitudes(4,iampl)) < tol8) then
1596 : 0 : sigma = amplitudes(5,iampl)
1597 : : end if
1598 : : end do
1599 : : !If not, the amplitude is zero
1600 [ # # # # : 0 : if(abs(sigma) < tol8)then
# # # ]
1601 : 0 : write (msg, '(a,I0,a,3es12.5,2a,I0)') ' The amplitude of the unstable mode ',&
1602 [ # # ]: 0 : int(imode),' of the qpt ',thm_scells(iconfig)%qphon(:), ch10,&
1603 : 0 : 'is set to zero for the configuration ',iconfig
1604 : 0 : ABI_WARNING(msg)
1605 : : end if
1606 : : end select
1607 : : end if
1608 : :
1609 : : ! add displacement for this mode to supercell positions eq 5 of Zacharias
1610 : 0 : call RANDOM_NUMBER(rand)
1611 : 0 : rand = two * rand - one
1612 : :
1613 : : ! from TDEP documentation for gaussian distribution of displacements
1614 : : !call RANDOM_NUMBER(rand1)
1615 : : !call RANDOM_NUMBER(rand2)
1616 : : ! rand = sqrt(-two*rand1) * sin(twopi*rand2)
1617 : :
1618 : : ! if (rand > half) then
1619 : : ! rand = one
1620 : : ! else
1621 : : ! rand = -one
1622 : : ! end if
1623 : :
1624 : 0 : freeze_displ = rand * sigma
1625 : :
1626 : 0 : call thm_scells(iconfig)%freeze_displ(phdispl1(:,:,:), freeze_displ)
1627 : : end do !iconfig
1628 : : end do !imode
1629 : : end do !iq
1630 : :
1631 [ # # ]: 0 : ABI_FREE(phfrq_allq)
1632 [ # # ]: 0 : ABI_FREE(phdispl_allq)
1633 [ # # ]: 0 : ABI_FREE(phdispl1)
1634 [ # # ]: 0 : ABI_FREE(qibz)
1635 [ # # ]: 0 : ABI_FREE(qbz)
1636 [ # # ]: 0 : ABI_FREE(wtqibz)
1637 : :
1638 [ # # # # : 0 : end subroutine thermal_supercell_make
# # ]
1639 : : !!***
1640 : :
1641 : : !----------------------------------------------------------------------
1642 : :
1643 : : !!****f* m_phonons/thermal_supercell_free
1644 : : !! NAME
1645 : : !! thermal_supercell_free
1646 : : !!
1647 : : !! FUNCTION
1648 : : !! deallocate thermal array of supercells
1649 : : !!
1650 : : !! SOURCE
1651 : :
1652 : 1 : subroutine thermal_supercell_free(nscells, thm_scells)
1653 : :
1654 : : !Arguments ------------------------------------
1655 : : !scalars
1656 : : integer, intent(in) :: nscells
1657 : : type(supercell_type), allocatable, intent(inout) :: thm_scells(:)
1658 : :
1659 : : ! local
1660 : : integer :: icell
1661 : :
1662 [ + - ]: 1 : if (allocated(thm_scells)) then
1663 [ + + ]: 11 : do icell = 1, nscells
1664 : 11 : call thm_scells(icell)%free()
1665 : : end do
1666 : : end if
1667 : :
1668 : 1 : end subroutine thermal_supercell_free
1669 : : !!***
1670 : :
1671 : : !----------------------------------------------------------------------
1672 : :
1673 : : !!****f* m_phonons/zacharias_supercell_print
1674 : : !! NAME
1675 : : !! zacharias_supercell_print
1676 : : !!
1677 : : !! FUNCTION
1678 : : !! print files with thermal array of supercells
1679 : : !!
1680 : : !! INPUTS
1681 : : !!
1682 : : !! OUTPUT
1683 : : !!
1684 : : !! SOURCE
1685 : :
1686 : 1 : subroutine zacharias_supercell_print(fname, ntemper, tempermin, temperinc, thm_scells)
1687 : :
1688 : : !Arguments ------------------------------------
1689 : : !scalars
1690 : : integer, intent(in) :: ntemper
1691 : : real(dp), intent(in) :: tempermin
1692 : : real(dp), intent(in) :: temperinc
1693 : : type(supercell_type), intent(in) :: thm_scells(ntemper)
1694 : : character(len=fnlen), intent(in) :: fname
1695 : :
1696 : : ! local
1697 : : integer :: itemp
1698 : : character(len=80) :: title1, title2
1699 : : character(len=fnlen) :: filename
1700 : : real(dp) :: temper
1701 : : character(len=10) :: temper_str
1702 : :
1703 [ + + ]: 11 : do itemp = 1, ntemper
1704 : 10 : temper = dble(itemp-1)*temperinc+tempermin
1705 : 10 : write(temper_str,'(I8)') int(temper)
1706 [ + - ]: 10 : write(filename, '(3a)') trim(fname), "_T_", trim(adjustl(temper_str))
1707 : 10 : write(title1, '(3a)') "# Zacharias thermalized supercell at temperature T= ", trim(temper_str), " Kelvin"
1708 : 10 : title2 = "# generated with alternating thermal displacements of all phonons"
1709 : 11 : call thm_scells(itemp)%print_abinit(filename, title1, title2)
1710 : : end do
1711 : :
1712 : 1 : end subroutine zacharias_supercell_print
1713 : : !!***
1714 : :
1715 : : !!****f* m_phonons/thermal_supercell_print
1716 : : !! NAME
1717 : : !! thermal_supercell_print
1718 : : !!
1719 : : !! FUNCTION
1720 : : !! print files with thermalized array of random supercell configurations
1721 : : !!
1722 : : !! INPUTS
1723 : : !!
1724 : : !! OUTPUT
1725 : : !!
1726 : : !! SOURCE
1727 : :
1728 : 0 : subroutine thermal_supercell_print(fname, nconfig, temperature_K, thm_scells)
1729 : :
1730 : : !Arguments ------------------------------------
1731 : : !scalars
1732 : : integer, intent(in) :: nconfig
1733 : : type(supercell_type), intent(in) :: thm_scells(nconfig)
1734 : : character(len=fnlen), intent(in) :: fname
1735 : : real(dp), intent(in) :: temperature_K
1736 : :
1737 : : ! local
1738 : : integer :: iconfig,itemp
1739 : : character(len=80) :: title1, title2
1740 : : character(len=fnlen) :: filename
1741 : : character(len=10) :: config_str
1742 : :
1743 [ # # ]: 0 : do iconfig = 1, nconfig
1744 : 0 : write(config_str,'(I8)') iconfig
1745 [ # # ]: 0 : write(filename, '(3a)') trim(fname), "_cf_", trim(adjustl(config_str))
1746 : 0 : write(title1, '(a,I6,a)') "# thermalized supercell at temperature T= ", temperature_K, " Kelvin"
1747 : 0 : title2 = "# generated with random thermal displacements of all phonons"
1748 : 0 : call thm_scells(itemp)%print_abinit(filename, title1, title2)
1749 : : end do
1750 : :
1751 : 0 : end subroutine thermal_supercell_print
1752 : : !!***
1753 : :
1754 : : !----------------------------------------------------------------------
1755 : :
1756 : : !!****f* m_phonons/phdos_ncwrite
1757 : : !! NAME
1758 : : !! phdos_ncwrite
1759 : : !!
1760 : : !! FUNCTION
1761 : : !! Save the content of the object in a netcdf file.
1762 : : !!
1763 : : !! INPUTS
1764 : : !! ncid=NC file handle (open in the caller)
1765 : : !! phdos<phdos_t>=Container object
1766 : : !!
1767 : : !! OUTPUT
1768 : : !! Only writing
1769 : : !!
1770 : : !! NOTES
1771 : : !! Frequencies are in eV, DOS are in states/eV.
1772 : : !!
1773 : : !! SOURCE
1774 : :
1775 : 63 : subroutine phdos_ncwrite(phdos, ncid)
1776 : :
1777 : : !Arguments ------------------------------------
1778 : : !scalars
1779 : : class(phdos_t),intent(in) :: phdos
1780 : : integer,intent(in) :: ncid
1781 : :
1782 : : !Local variables-------------------------------
1783 : : integer :: ncerr
1784 : : ! *************************************************************************
1785 : :
1786 : : ! Define dimensions
1787 [ - + ]: 63 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
1788 : :
1789 : : ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_atoms", phdos%natom),&
1790 : : nctkdim_t("number_of_atom_species", phdos%ntypat), nctkdim_t("number_of_frequencies", phdos%nomega), &
1791 [ + + ]: 378 : nctkdim_t("nqibz", phdos%nqibz), nctkdim_t("number_of_normal_vec_dmm", phdos%n_normal_vec_dmm)])
1792 [ - + ]: 63 : NCF_CHECK(ncerr)
1793 : :
1794 : : !scalars
1795 [ + + - + : 126 : NCF_CHECK(nctk_def_iscalars(ncid, ["prtdos"]))
- - ]
1796 [ + + - + : 126 : NCF_CHECK(nctk_def_dpscalars(ncid, ["dossmear"]))
- - ]
1797 : :
1798 : : !arrays
1799 : : ncerr = nctk_def_arrays(ncid, [&
1800 : : nctkarr_t('wmesh', "dp", 'number_of_frequencies'),&
1801 : : nctkarr_t('phdos', "dp", 'number_of_frequencies'),&
1802 : : nctkarr_t('phdos_dmm', "dp", 'number_of_frequencies, number_of_normal_vec_dmm'),&
1803 : : nctkarr_t('normal_vec_dmm', "dp", 'three, number_of_normal_vec_dmm'),&
1804 : : nctkarr_t('pjdos', "dp", 'number_of_frequencies, three, number_of_atoms'),&
1805 : : nctkarr_t('pjdos_type', "dp", 'number_of_frequencies, number_of_atom_species'),&
1806 : : nctkarr_t('pjdos_rc_type', "dp", 'number_of_frequencies, three, number_of_atom_species'), &
1807 : : nctkarr_t('msqd_dos_atom', "dp", 'number_of_frequencies, three, three, number_of_atoms'), &
1808 : : nctkarr_t('qptrlatt', "int", 'three, three'), &
1809 : : nctkarr_t('shiftq', "dp", 'three') &
1810 [ + + ]: 693 : ])
1811 [ - + ]: 63 : NCF_CHECK(ncerr)
1812 : :
1813 : : ! Write variables. Note unit conversion.
1814 [ - + ]: 63 : NCF_CHECK(nctk_set_datamode(ncid))
1815 [ - + ]: 63 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "prtdos"), phdos%prtdos))
1816 [ - + ]: 63 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'dossmear'), phdos%dossmear*Ha_eV))
1817 [ + - + + : 85905 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'wmesh'), phdos%omega*Ha_eV))
- + - - -
- ]
1818 [ - + ]: 63 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'normal_vec_dmm'), phdos%normal_vec_dmm))
1819 [ + - + + : 85905 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdos'), phdos%phdos/Ha_eV))
- + - - -
- ]
1820 [ + - + + : 601398 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdos_dmm'), phdos%phdos_dmm/Ha_eV)) ! TODO: these should also be converted for the velocity
+ + - + -
- - - -
- ]
1821 [ + - + + : 592900 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos'), phdos%pjdos/Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1822 [ + - + + : 122611 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos_type'), phdos%pjdos_type/Ha_eV))
+ + - + -
- - - -
- ]
1823 [ + - + + : 367802 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'pjdos_rc_type'), phdos%pjdos_rc_type/Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1824 [ + - + + : 1778695 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'msqd_dos_atom'), phdos%msqd_dos_atom/Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1825 [ - + ]: 63 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qptrlatt'), phdos%qptrlatt))
1826 [ - + ]: 63 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'shiftq'), phdos%shiftq))
1827 : :
1828 : 63 : end subroutine phdos_ncwrite
1829 : : !!***
1830 : :
1831 : : !----------------------------------------------------------------------
1832 : :
1833 : : !!****f* m_phonons/mkphbs
1834 : : !! NAME
1835 : : !! mkphbs
1836 : : !!
1837 : : !! FUNCTION
1838 : : !! Function to calculate the phonon band structure, from the IFC
1839 : : !!
1840 : : !! INPUTS
1841 : : !! Ifc<ifc_type>=Interatomic force constants
1842 : : !! crystal<type(crystal_t)> = Info on the crystalline structure.
1843 : : !! inp= (derived datatype) contains all the input variables
1844 : : !! ddb<type(ddb_type)>=Object storing the DDB results.
1845 : : !! asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
1846 : : !! prefix=Prefix for output files.
1847 : : !! dielt(3,3)=dielectric tensor
1848 : : !! comm=MPI communicator
1849 : : !!
1850 : : !! OUTPUT
1851 : : !! Only writing.
1852 : : !!
1853 : : !! SOURCE
1854 : :
1855 : 65 : subroutine mkphbs(Ifc,Crystal,inp,ddb,asrq0,prefix,comm)
1856 : :
1857 : : !Arguments -------------------------------
1858 : : !scalars
1859 : : integer,intent(in) :: comm
1860 : : character(len=*),intent(in) :: prefix
1861 : : type(ifc_type),intent(in) :: Ifc
1862 : : type(crystal_t),intent(in) :: Crystal
1863 : : type(anaddb_dataset_type),target,intent(in) :: inp
1864 : : type(ddb_type),intent(in) :: ddb
1865 : : type(asrq0_t),intent(inout) :: asrq0
1866 : : !Local variables -------------------------
1867 : : !scalars
1868 : : integer,parameter :: master=0
1869 : : integer :: unt, iphl1,iblok,rftyp, ii,nfineqpath,nsym,natom,ncid,nprocs,my_rank
1870 : : integer :: natprj_bs,eivec,enunit,ifcflag,ptgroupma,spgroup
1871 : : real(dp) :: freeze_displ, cfact, omega, omega_min, gaussmaxarg, gaussfactor, gaussprefactor, xx, eta
1872 : : character(500) :: msg
1873 : : character(len=8) :: unitname
1874 : : !arrays
1875 : : integer :: bravais(11),rfphon(4),rfelfd(4),rfstrs(4), units(2)
1876 : : integer :: nomega, imode, iomega
1877 : 65 : integer,allocatable :: ndiv(:)
1878 : : real(dp) :: speedofsound(3),genafm(3)
1879 : : real(dp) :: qphnrm(3), qphon(3), qphon_padded(3,3),res(3)
1880 : 130 : real(dp) :: d2cart(2,ddb%msize),real_qphon(3)
1881 : 130 : real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom),eigval(3,Crystal%natom),phangmom(3,3*Crystal%natom)
1882 : 65 : real(dp),allocatable :: phfrq(:),eigvec(:,:,:,:,:)
1883 : 65 : real(dp),allocatable :: save_phfrq(:,:),save_phdispl_cart(:,:,:,:),save_qpoints(:,:),save_phangmom(:,:,:)
1884 : 65 : real(dp),allocatable :: weights(:), dos4bs(:)
1885 : 65 : real(dp),allocatable,target :: alloc_path(:,:)
1886 : 65 : real(dp),pointer :: fineqpath(:,:)
1887 : 65 : type(atprj_type) :: atprj
1888 : : ! *********************************************************************
1889 : :
1890 : : ! Only master works for the time being
1891 : 65 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1892 [ - + ]: 65 : if (my_rank /= master) return
1893 : :
1894 [ + + ]: 195 : units = [std_out, ab_out]
1895 : 65 : nsym = Crystal%nsym; natom = Crystal%natom
1896 : :
1897 : : ! Copy parameters from inp (then I will try to remove inp from the API so that I can call mkphbs in eph)
1898 : 65 : ifcflag = inp%ifcflag
1899 : 65 : natprj_bs = inp%natprj_bs
1900 : 65 : freeze_displ = inp%freeze_displ
1901 : 65 : eivec = inp%eivec; enunit = inp%enunit
1902 : 65 : rftyp=inp%rfmeth
1903 : :
1904 : 65 : nullify(fineqpath)
1905 : 65 : nfineqpath = inp%nph1l
1906 : 65 : fineqpath => inp%qph1l
1907 : :
1908 [ + + ]: 65 : if(inp%nph1l==0) then
1909 [ - + ]: 8 : if (inp%nqpath==0) then
1910 : : return ! if there is nothing to do, return
1911 : : else
1912 : : ! allow override of nph1l with nqpath if the former is not set
1913 : : ! allocate and compute path here and make fineqpath points to it
1914 [ + - - + : 24 : ABI_MALLOC(ndiv,(inp%nqpath-1))
+ - + - +
- - + -
+ ]
1915 : 8 : call make_path(inp%nqpath,inp%qpath,Crystal%gmet,'G',inp%ndivsm,ndiv,nfineqpath,alloc_path,std_out)
1916 [ - + ]: 8 : ABI_FREE(ndiv)
1917 : 8 : fineqpath => alloc_path
1918 : : end if
1919 : : end if
1920 : :
1921 [ + - + + ]: 5265 : write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10
1922 : 65 : call wrtout(units, msg)
1923 : :
1924 [ + + ]: 65 : if (natprj_bs > 0) call atprj%init(natom, natprj_bs, inp%iatprj_bs, prefix)
1925 : :
1926 [ + - - + : 260 : ABI_MALLOC(phfrq, (3*natom))
+ - + - +
- - + -
+ ]
1927 [ - + + - : 325 : ABI_MALLOC(eigvec, (2,3,natom,3,natom))
+ - - + +
- + - + -
- + - + ]
1928 [ - + + - : 195 : ABI_MALLOC(save_qpoints, (3,nfineqpath))
+ - + - -
+ - + ]
1929 [ - + + - : 260 : ABI_MALLOC(save_phfrq, (3*natom,nfineqpath))
- + + - +
- + - - +
- + ]
1930 [ - + + - : 325 : ABI_MALLOC(save_phdispl_cart, (2,3*natom,3*natom,nfineqpath))
- + + - -
+ + - + -
+ - - + -
+ ]
1931 [ - + + - : 260 : ABI_MALLOC(save_phangmom, (3,3*natom,nfineqpath))
- + + - +
- + - - +
- + ]
1932 [ + + ]: 260 : qphnrm = one
1933 : :
1934 [ + + ]: 1309 : do iphl1=1,nfineqpath
1935 : :
1936 : : ! Initialisation of the phonon wavevector
1937 [ + + ]: 4976 : qphon(:)=fineqpath(:,iphl1)
1938 : :
1939 [ + + ]: 1244 : if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1)
1940 : :
1941 [ + + ]: 4976 : save_qpoints(:,iphl1) = qphon / qphnrm(1)
1942 : :
1943 : : ! Generation of the dynamical matrix in cartesian coordinates
1944 [ + + ]: 1244 : if (ifcflag == 1) then
1945 : :
1946 : : ! Get phonon frequencies and displacements in reduced coordinates for this q-point
1947 : : !call ifc%fourq(cryst, save_qpoints(:,iphl1), phfrq, displ, out_eigvec=eigvec)
1948 : :
1949 : : ! Get d2cart using the interatomic forces and the
1950 : : ! long-range coulomb interaction through Ewald summation
1951 : : call gtdyn9(Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,Ifc%gprimd,Ifc%mpert,natom, &
1952 : : Ifc%nrpt,qphnrm(1),qphon,Crystal%rmet,Ifc%rprimd,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,ifc%zeff,&
1953 : : ifc%qdrp_cart,ifc%ewald_option,eta,xmpi_comm_self,Ifc%sys_dim,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad,&
1954 : 1202 : dielt_env=Ifc%dielt_env,dielt_thick=Ifc%dielt_thick)
1955 [ + + ]: 1202 : if (asrq0%asr==6) then
1956 [ + + ]: 644 : qphon_padded = zero; qphon_padded(:,1) = qphon(:)
1957 : 161 : call asrq0%apply(natom, ddb%mpert, ddb%msize, qphon_padded, Crystal, d2cart)
1958 : : end if
1959 [ + - ]: 42 : else if (ifcflag == 0) then
1960 : :
1961 : : !call ddb_diagoq(ddb, crystal, save_qpoints(:,iphl1), asrq0, ifc%symdynmat, rftyp, phfrq, displ, &
1962 : : ! out_eigvec=eigvec)
1963 : :
1964 : : ! Look for the information in the DDB (no interpolation here!)
1965 [ + + + + : 294 : rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0
+ + ]
1966 [ + + ]: 168 : qphon_padded = zero; qphon_padded(:,1) = qphon
1967 : :
1968 : 42 : call ddb%get_block(iblok,qphon_padded,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
1969 : :
1970 : : ! Copy the dynamical matrix in d2cart
1971 [ + + + + ]: 1159908 : d2cart(:,1:ddb%msize)=ddb%val(:,:,iblok)
1972 : :
1973 : : ! Eventually impose the acoustic sum rule based on previously calculated d2asr
1974 : 42 : call asrq0%apply(natom, ddb%mpert, ddb%msize, qphon_padded, Crystal, d2cart)
1975 : : end if
1976 : :
1977 : : ! Use inp%symdynmat instead of ifc because of ifcflag
1978 : : ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
1979 : : call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,&
1980 : : ddb%mpert,Crystal%nsym,natom,nsym,Crystal%ntypat,phfrq,qphnrm(1),qphon,&
1981 : 1244 : crystal%rprimd,inp%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
1982 : :
1983 : : ! Calculation of the phonon angular momentum
1984 : : ! maybe add it in dpft_phfrq directly ?
1985 : 1244 : call phangmom_from_eigvec(natom, eigvec, phangmom)
1986 : :
1987 [ + + ]: 1244 : if (abs(freeze_displ) > tol10) then
1988 : 6 : real_qphon = zero
1989 [ + - + + ]: 24 : if (abs(qphnrm(1)) > tol8) real_qphon = qphon / qphnrm(1)
1990 : : call freeze_displ_allmodes(displ, freeze_displ, natom, prefix, phfrq, &
1991 : 6 : real_qphon, crystal%rprimd, Crystal%typat, crystal%xcart, crystal%znucl)
1992 : : end if
1993 : :
1994 : : ! If requested, output projection of each mode on given atoms
1995 [ + + ]: 1244 : if (natprj_bs > 0) call atprj%print(iphl1, phfrq, eigvec)
1996 : :
1997 : : ! In case eivec == 4, write output files for band2eps (visualization of phonon band structures)
1998 [ + + ]: 1244 : if (eivec == 4) call sortph(eigvec,displ,strcat(prefix, "_B2EPS"),natom,phfrq)
1999 : :
2000 : : ! Write the phonon frequencies
2001 : 1244 : call dfpt_prtph(displ,eivec,enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
2002 : :
2003 [ + + ]: 14621 : save_phfrq(:,iphl1) = phfrq
2004 [ + + ]: 4976 : save_phdispl_cart(:,:,:,iphl1) = RESHAPE(displ, [2, 3*natom, 3*natom])
2005 [ + + ]: 3732 : save_phangmom(:,:,iphl1) = RESHAPE(phangmom, [3, 3*natom])
2006 : :
2007 : : ! Determine the symmetries of the phonon mode at Gamma
2008 : : ! TODO: generalize for other q-point little groups.
2009 [ + + + + ]: 4976 : if (sum(abs(qphon)) < DDB_QTOL) then
2010 : : call symanal(bravais,0,genafm,nsym,nsym,ptgroupma,Crystal%rprimd,spgroup, &
2011 : 72 : Crystal%symafm,Crystal%symrel,Crystal%tnons,tol5,verbose=.TRUE.)
2012 : 72 : call dfpt_symph(ab_out,ddb%acell,eigvec,Crystal%indsym,natom,nsym,phfrq,ddb%rprim,Crystal%symrel)
2013 : : end if
2014 : :
2015 : : ! if we have an acoustic mode (small q and acoustic type displacements)
2016 : : ! extrapolate speed of sound in this direction, and Debye frequency
2017 [ + + ]: 4976 : call wrap2_pmhalf(qphon, real_qphon, res)
2018 [ + + ]: 311 : if (sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) < quarter .and. &
2019 [ + + ]: 1309 : sqrt(real_qphon(1)**2+real_qphon(2)**2+real_qphon(3)**2) > tol6) then
2020 : 237 : call phdos_calc_vsound(eigvec, Crystal%gmet, natom, phfrq, real_qphon, speedofsound)
2021 : 237 : if (my_rank == master) call phdos_print_vsound(ab_out, Crystal%ucvol, speedofsound)
2022 : : end if
2023 : :
2024 : : end do ! iphl1
2025 : :
2026 : : ! calculate dos for the specific q points along the BS calculated. only Gaussians are possible - no interpolation
2027 [ + + + + : 14686 : omega_min = minval(save_phfrq(:,:))
+ + + - -
+ - - ]
2028 [ + + + + : 14686 : nomega=NINT( (maxval(save_phfrq(:,:))-omega_min) / inp%dosdeltae ) + 1
+ + + - -
+ - - ]
2029 : 65 : nomega=MAX(6,nomega) ! Ensure Simpson integration will be ok
2030 : :
2031 [ - + + - : 249018 : ABI_CALLOC(dos4bs,(nomega))
+ - - + -
+ + + ]
2032 : :
2033 : 65 : gaussmaxarg = sqrt(-log(1.d-90))
2034 : 65 : gaussprefactor = one/(inp%dossmear*sqrt(two_pi))
2035 : 65 : gaussfactor = one/(sqrt2*inp%dossmear)
2036 [ + + ]: 1309 : do iphl1=1,nfineqpath
2037 [ + + ]: 14686 : do imode=1,3*natom
2038 [ + + ]: 55778618 : do iomega=1, nomega
2039 : 55763997 : omega = omega_min + (iomega-1) * inp%dosdeltae
2040 : 55763997 : xx = (omega - save_phfrq(imode,iphl1)) * gaussfactor
2041 [ + + ]: 55777374 : if(abs(xx) < gaussmaxarg) dos4bs(iomega) = dos4bs(iomega) + gaussprefactor*exp(-xx*xx)
2042 : : end do
2043 : : end do
2044 : : end do
2045 : :
2046 : : !deallocate sortph array
2047 : 65 : call end_sortph()
2048 : :
2049 [ + + ]: 65 : if (natprj_bs > 0) call atprj%free()
2050 : :
2051 : : ! WRITE OUT FILES
2052 : : if (my_rank == master) then
2053 [ - + + - : 195 : ABI_MALLOC(weights, (nfineqpath))
+ - + - -
+ - + ]
2054 [ + + ]: 1309 : weights = one
2055 [ - + ]: 65 : NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST")
2056 [ - + ]: 65 : NCF_CHECK(crystal%ncwrite(ncid))
2057 : 65 : call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart,save_phangmom)
2058 : : ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
2059 [ + + + + ]: 65 : if (inp%nph2l /= 0 .and. inp%ifcflag == 1) then
2060 : 14 : call ifc%calcnwrite_nana_terms(crystal, inp%nph2l, inp%qph2l, inp%qnrml2, ncid)
2061 : : end if
2062 [ - + ]: 65 : NCF_CHECK(nf90_close(ncid))
2063 : :
2064 : 65 : call phonons_write_phfrq(prefix, natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart, save_phangmom)
2065 : :
2066 : 65 : select case (inp%prtphbands)
2067 : : case (0)
2068 : 65 : continue
2069 : :
2070 : : case (1)
2071 [ + + ]: 65 : if (inp%nph1l == 0) then
2072 : : call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq, &
2073 : 8 : qptbounds=inp%qpath)
2074 : : else
2075 : 57 : call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nfineqpath, save_qpoints, save_phfrq)
2076 : : end if
2077 : :
2078 : : case (2)
2079 [ # # ]: 0 : if (inp%nph1l == 0) then
2080 : 0 : call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq, qptbounds=inp%qpath)
2081 : : else
2082 : 0 : call phonons_write_gnuplot(prefix, natom, nfineqpath, save_qpoints, save_phfrq)
2083 : : end if
2084 : :
2085 : : case default
2086 [ + - - - ]: 65 : ABI_WARNING(sjoin("Don't know how to handle prtphbands:", itoa(inp%prtphbands)))
2087 : : end select
2088 : :
2089 : : ! write out DOS file for q along this path
2090 : 65 : cfact=one
2091 : 65 : unitname = 'Ha'
2092 [ - + ]: 65 : if (open_file('PHBST_partial_DOS',msg,newunit=unt,form="formatted",action="write") /= 0) then
2093 : 0 : ABI_ERROR(msg)
2094 : : end if
2095 : 65 : write(msg,'(3a)')'# ',ch10,'# Partial phonon density of states for q along a band structure path'
2096 : 65 : call wrtout(unt, msg)
2097 : 65 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
2098 : 65 : call wrtout(unt, msg)
2099 : 65 : write(msg,'(a,E20.10,2a,i8)') '# Gaussian method with smearing = ',inp%dossmear*cfact,unitname, ', nq =', nfineqpath
2100 : 65 : call wrtout(unt, msg)
2101 : 65 : write(msg,'(5a)')'# ',ch10,'# omega PHDOS ',ch10,'# '
2102 : 65 : call wrtout(unt, msg)
2103 [ + + ]: 248888 : do iomega=1,nomega
2104 : 248823 : omega = omega_min + (iomega-1) * inp%dosdeltae
2105 : 248823 : write(unt,'(2es17.8)',advance='NO')omega*cfact,dos4bs(iomega)/cfact
2106 : 248888 : write(unt,*)
2107 : : end do
2108 : 65 : close(unt)
2109 : :
2110 [ - + ]: 65 : ABI_FREE(weights)
2111 : : end if
2112 : :
2113 [ - + ]: 65 : ABI_FREE(save_qpoints)
2114 [ - + ]: 65 : ABI_FREE(save_phfrq)
2115 [ - + ]: 65 : ABI_FREE(save_phdispl_cart)
2116 [ - + ]: 65 : ABI_FREE(save_phangmom)
2117 [ - + ]: 65 : ABI_FREE(phfrq)
2118 [ - + ]: 65 : ABI_FREE(eigvec)
2119 [ - + ]: 65 : ABI_FREE(dos4bs)
2120 [ + + ]: 65 : ABI_SFREE(alloc_path)
2121 : :
2122 [ - + - + : 130 : end subroutine mkphbs
- + ]
2123 : : !!***
2124 : :
2125 : : !!****f* m_phonons/phdos_calc_vsound
2126 : : !!
2127 : : !! NAME
2128 : : !! phdos_calc_vsound
2129 : : !!
2130 : : !! FUNCTION
2131 : : !! From the frequencies for acoustic modes at small q, estimate speed of sound (which also gives Debye temperature)
2132 : : !!
2133 : : !! INPUTS
2134 : : !! eigvec(2,3*natom,3*natom) = phonon eigenvectors at present q-point
2135 : : !! gmet(3,3) = metric tensor in reciprocal space.
2136 : : !! natom = number of atoms in the unit cell
2137 : : !! phfrq(3*natom) = phonon frequencies at present q-point
2138 : : !! qphon(3) = phonon q-point
2139 : : !! ucvol = unit cell volume
2140 : : !!
2141 : : !! OUTPUT
2142 : : !!
2143 : : !! SOURCE
2144 : :
2145 : 10606 : subroutine phdos_calc_vsound(eigvec, gmet, natom, phfrq, qphon, speedofsound)
2146 : :
2147 : : !Arguments -------------------------------
2148 : : !scalars
2149 : : integer, intent(in) :: natom
2150 : : !arrays
2151 : : real(dp), intent(in) :: gmet(3,3),qphon(3),phfrq(3*natom),eigvec(2,3*natom,3*natom)
2152 : : real(dp), intent(out) :: speedofsound(3)
2153 : :
2154 : : !Local variables -------------------------
2155 : : integer :: iatref,imode, iatom, isacoustic, imode_acoustic
2156 : : ! character(len=500) :: msg
2157 : : real(dp) :: qnormcart
2158 : : real(dp) :: qtmp(3)
2159 : : ! *********************************************************************
2160 : :
2161 : 10606 : imode_acoustic = 0
2162 : :
2163 [ + + ]: 66973 : do imode = 1, 3*natom
2164 : : ! Check if this mode is acoustic like: scalar product of all displacement vectors are collinear
2165 : 179172 : isacoustic = 1
2166 : : ! Find reference atom with non-zero displacement
2167 [ + + ]: 179172 : do iatom=1,natom
2168 [ + + + + : 1284417 : if(sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3,imode)**2) >tol16)iatref=iatom
+ + ]
2169 : : enddo
2170 : : ! Now compute scalar product, and check they are all positive
2171 [ + + ]: 179172 : do iatom = 1, natom
2172 [ + + + + : 1228050 : if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, imode)&
+ + ]
2173 : 88900 : *eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, imode)) < tol16 ) isacoustic = 0
2174 : : end do
2175 [ + + ]: 56367 : if (isacoustic == 0) cycle
2176 : 31797 : imode_acoustic = min(imode_acoustic + 1, 3)
2177 : :
2178 : : ! write (msg, '(a,I6,a,3F12.4)') ' Found acoustic mode ', imode, ' for |q| in red coord < 0.25 ; q = ', qphon
2179 : : ! call wrtout(std_out, msg)
2180 [ + + + + ]: 413361 : qtmp = matmul(gmet, qphon)
2181 [ + + ]: 127188 : qnormcart = two * pi * sqrt(sum(qphon*qtmp))
2182 : 66973 : speedofsound(imode_acoustic) = phfrq(imode) / qnormcart
2183 : : end do
2184 : :
2185 : 10606 : end subroutine phdos_calc_vsound
2186 : : !!***
2187 : :
2188 : : !!****f* m_phonons/phdos_print_vsound
2189 : : !!
2190 : : !! NAME
2191 : : !! phdos_print_vsound
2192 : : !!
2193 : : !! FUNCTION
2194 : : !! Print out estimate speed of sound and Debye temperature at this (small) q
2195 : : !! should only be called by master proc for the hard unit number
2196 : : !!
2197 : : !! INPUTS
2198 : : !! unit=Fortran unit number
2199 : : !! speedofsound(3)
2200 : : !!
2201 : : !! OUTPUT
2202 : : !! Only writing
2203 : : !!
2204 : : !! SOURCE
2205 : :
2206 : 237 : subroutine phdos_print_vsound(iunit, ucvol, speedofsound)
2207 : :
2208 : : !Arguments -------------------------------
2209 : : !scalras
2210 : : integer, intent(in) :: iunit
2211 : : real(dp), intent(in) :: ucvol
2212 : : !arrays
2213 : : real(dp), intent(in) :: speedofsound(3)
2214 : :
2215 : : !Local variables -------------------------
2216 : : integer :: imode_acoustic, units(2)
2217 : : character(len=500) :: msg
2218 : : real(dp) :: tdebye
2219 : : ! *********************************************************************
2220 : :
2221 [ + + ]: 711 : units = [std_out, iunit]
2222 : :
2223 [ + + ]: 948 : do imode_acoustic = 1, 3
2224 : : ! from phonon frequency, estimate speed of sound by linear interpolation from Gamma
2225 : : write (msg, '(2a,a,E20.10,a,a,F20.5)') &
2226 : 711 : ' Speed of sound for this q and mode:',ch10,&
2227 : 711 : ' in atomic units: ', speedofsound(imode_acoustic), ch10,&
2228 : 1422 : ' in units km/s: ', speedofsound(imode_acoustic) * Bohr_Ang * 1.d-13 / Time_Sec
2229 : 711 : call wrtout(units, msg)
2230 : :
2231 : : ! also estimate partial Debye temperature, = energy if this band went to zone edge
2232 : 711 : tdebye = speedofsound(imode_acoustic) * pi * (six / pi / ucvol)**(third)
2233 : : write (msg, '(2a,a,E20.10,a,a,F20.5)') &
2234 : 711 : ' Partial Debye temperature for this q and mode:',ch10,&
2235 : 711 : ' in atomic units: ', tdebye, ch10,&
2236 : 1422 : ' in SI units K : ', tdebye * Ha_K
2237 : 711 : call wrtout(units, msg)
2238 : 948 : call wrtout(units, "")
2239 : : end do
2240 : :
2241 : 237 : end subroutine phdos_print_vsound
2242 : : !!***
2243 : :
2244 : : !----------------------------------------------------------------------
2245 : :
2246 : : !!****f* m_phonons/phdos_print_msqd
2247 : : !!
2248 : : !! NAME
2249 : : !! phdos_print_msqd
2250 : : !!
2251 : : !! FUNCTION
2252 : : !! Print out mean square displacement and velocity for each atom (trace and full matrix) as a function of T
2253 : : !! see for example https://atztogo.github.io/phonopy/thermal-displacement.html#thermal-displacement
2254 : : !! Only master node should call this routine.
2255 : : !!
2256 : : !! INPUTS
2257 : : !! PHdos structure
2258 : : !!
2259 : : !! OUTPUT
2260 : : !! to file only
2261 : : !!
2262 : : !! SOURCE
2263 : :
2264 : 11 : subroutine phdos_print_msqd(PHdos, fname, ntemper, tempermin, temperinc)
2265 : :
2266 : : !Arguments -------------------------------
2267 : : !scalars
2268 : : class(phdos_t),intent(in) :: PHdos
2269 : : character(len=*),intent(in) :: fname
2270 : : integer, intent(in) :: ntemper
2271 : : real(dp), intent(in) :: tempermin, temperinc
2272 : :
2273 : : !Local variables -------------------------
2274 : : integer :: io, iomin, itemp, iunit, junit, iatom
2275 : : real(dp) :: temper
2276 : : character(len=500) :: msg
2277 : : character(len=fnlen) :: fname_msqd, fname_veloc
2278 : : !arrays
2279 : 11 : real(dp), allocatable :: bose_msqd(:,:), tmp_msqd(:,:), integ_msqd(:,:)
2280 : 11 : real(dp), allocatable :: bose_msqv(:,:), tmp_msqv(:,:), integ_msqv(:,:)
2281 : : ! *********************************************************************
2282 : :
2283 [ + - ]: 11 : fname_msqd = trim(fname) //"_MSQD_T"
2284 [ - + ]: 11 : if (open_file(fname_msqd, msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2285 : 0 : ABI_ERROR(msg)
2286 : : end if
2287 [ + - ]: 11 : fname_veloc = trim(fname) // "_MSQV_T"
2288 [ - + ]: 11 : if (open_file(fname_veloc, msg, newunit=junit, form="formatted", status="unknown", action="write") /= 0) then
2289 : 0 : ABI_ERROR(msg)
2290 : : end if
2291 : :
2292 : : ! write the header
2293 : 11 : write (iunit, '(a)') '# mean square displacement for each atom as a function of T (bohr^2)'
2294 : 11 : write (junit, '(a)') "# mean square velocity for each atom as a function of T (bohr^2/atomic time unit^2)"
2295 : :
2296 : 11 : write (msg, '(a,F18.10,a,F18.10,a)') '# T in Kelvin, from ', tempermin, ' to ', tempermin+(ntemper-1)*temperinc
2297 : 11 : write (iunit, '(a)') trim(msg)
2298 : 11 : write (junit, '(a)') trim(msg)
2299 : :
2300 : 11 : write (msg, '(2a)') '# T |u^2| u_xx u_yy u_zz',&
2301 : 22 : ' u_yz u_xz u_xy in bohr^2'
2302 : 11 : write (iunit, '(a)') trim(msg)
2303 : 11 : write (msg, '(3a)') '# T |v^2| v_xx v_yy v_zz',&
2304 : 11 : ' v_yz v_xz v_xy',&
2305 : 22 : ' in bohr^2/atomic time unit^2'
2306 : 11 : write (junit, '(a)') trim(msg)
2307 : :
2308 [ - + + - : 33 : ABI_MALLOC(tmp_msqd, (PHdos%nomega,9))
+ - + - -
+ - + ]
2309 [ - + + - : 22 : ABI_MALLOC(tmp_msqv, (PHdos%nomega,9))
- + - + ]
2310 [ - + + - : 33 : ABI_MALLOC(integ_msqd, (9,ntemper))
+ - + - -
+ - + ]
2311 [ - + + - : 22 : ABI_MALLOC(integ_msqv, (9,ntemper))
- + - + ]
2312 [ - + + - : 44 : ABI_MALLOC(bose_msqd, (PHdos%nomega, ntemper))
- + + - +
- + - - +
- + ]
2313 [ - + + - : 33 : ABI_MALLOC(bose_msqv, (PHdos%nomega, ntemper))
- + + - -
+ - + ]
2314 : :
2315 [ - + ]: 10455 : do io=1, PHdos%nomega
2316 [ + + ]: 10455 : if ( PHdos%omega(io) >= 2._dp * 4.56d-6 ) exit ! 2 cm-1 TODO: make this an input parameter
2317 : : end do
2318 : : iomin = io
2319 : :
2320 : : ! calculate bose only once for each atom (instead of for each atom)
2321 [ + + + + ]: 407961 : bose_msqd = zero
2322 [ + + + + ]: 407961 : bose_msqv = zero
2323 [ + + ]: 121 : do itemp = 1, ntemper
2324 : 110 : temper = tempermin + (itemp-1) * temperinc
2325 [ + - ]: 110 : if (temper < 1.e-3) cycle ! millikelvin at least to avoid exploding Bose factor(TM)
2326 [ + + ]: 303521 : do io = iomin, PHdos%nomega
2327 : : ! NB: factors follow convention in phonopy documentation
2328 : : ! the 1/sqrt(omega) factor in phonopy is contained in the displacement vector definition
2329 : : ! bose() is dimensionless
2330 : : !bose_msqd(io, itemp) = (half + one / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) / PHdos%omega(io)
2331 : : !bose_msqv(io, itemp) = (half + one / ( exp(PHdos%omega(io)/(kb_HaK*temper)) - one )) * PHdos%omega(io)
2332 : 303400 : bose_msqd(io, itemp) = (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) / PHdos%omega(io)
2333 : 303510 : bose_msqv(io, itemp) = (half + bose_einstein(PHdos%omega(io),kb_HaK*temper)) * PHdos%omega(io)
2334 : : end do
2335 : : end do
2336 : :
2337 [ + + ]: 36 : do iatom=1,PHdos%natom
2338 : 25 : write (msg, '(a,I8)') '# atom number ', iatom
2339 : 25 : write (iunit, '(a)') trim(msg)
2340 : 25 : write (junit, '(a)') trim(msg)
2341 : :
2342 : : ! for each T and each atom, integrate msqd matrix with Bose Einstein factor and output
2343 [ + + + + ]: 2525 : integ_msqd = zero
2344 [ + + ]: 75 : tmp_msqd = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/))
2345 : :
2346 : : ! perform all integrations as matrix multiplication: integ_msqd (idir, itemp) = [tmp_msqd(io,idir)]^T * bose_msqd(io,itemp)
2347 : 25 : call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqd,PHdos%nomega, bose_msqd, PHdos%nomega, zero, integ_msqd, 9)
2348 : : ! NB: this presumes an equidistant omega grid
2349 [ + + + + ]: 2525 : integ_msqd = integ_msqd * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom)
2350 : :
2351 [ + + + + ]: 2525 : integ_msqv = zero
2352 [ + + ]: 75 : tmp_msqv = reshape(PHdos%msqd_dos_atom(:,:,:,iatom), (/PHdos%nomega, 9/))
2353 : :
2354 : : ! perform all integrations as matrix multiplication: integ_msqv (idir, itemp) = [tmp_msqv(io,idir)]^T * bose_msqv(io,itemp)
2355 : 25 : call DGEMM('T','N', 9, ntemper, PHdos%nomega, one, tmp_msqv,PHdos%nomega, bose_msqv, PHdos%nomega, zero, integ_msqv, 9)
2356 : : ! NB: this presumes an equidistant omega grid
2357 [ + + + + ]: 2525 : integ_msqv = integ_msqv * (PHdos%omega(2)-PHdos%omega(1)) / PHdos%atom_mass(iatom)
2358 : :
2359 : : ! print out stuff
2360 [ + + ]: 275 : do itemp = 1, ntemper
2361 : 250 : temper = tempermin + (itemp-1) * temperinc
2362 : : write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') &
2363 : 250 : temper, third*(integ_msqd(1,itemp)+integ_msqd(5,itemp)+integ_msqd(9,itemp)), &
2364 : 250 : integ_msqd(1,itemp),integ_msqd(5,itemp),integ_msqd(9,itemp), &
2365 : 500 : integ_msqd(6,itemp),integ_msqd(3,itemp),integ_msqd(2,itemp)
2366 : 250 : write (iunit, '(a)') trim(msg)
2367 : : write (msg, '(F10.2,4x,E22.10,2x,6E22.10)') &
2368 : 250 : temper, third*(integ_msqv(1,itemp)+integ_msqv(5,itemp)+integ_msqv(9,itemp)), &
2369 : 250 : integ_msqv(1,itemp),integ_msqv(5,itemp),integ_msqv(9,itemp), &
2370 : 500 : integ_msqv(6,itemp),integ_msqv(3,itemp),integ_msqv(2,itemp)
2371 : 275 : write (junit, '(a)') trim(msg)
2372 : : end do ! itemp
2373 : :
2374 : 25 : write (iunit, '(a)') ''
2375 : 36 : write (junit, '(a)') ''
2376 : : enddo ! iatom
2377 : :
2378 [ - + ]: 11 : ABI_FREE(tmp_msqd)
2379 [ - + ]: 11 : ABI_FREE(tmp_msqv)
2380 [ - + ]: 11 : ABI_FREE(bose_msqd)
2381 [ - + ]: 11 : ABI_FREE(bose_msqv)
2382 [ - + ]: 11 : ABI_FREE(integ_msqd)
2383 [ - + ]: 11 : ABI_FREE(integ_msqv)
2384 : :
2385 : 11 : close(iunit)
2386 : 11 : close(junit)
2387 : :
2388 : 11 : end subroutine phdos_print_msqd
2389 : : !!***
2390 : :
2391 : : !----------------------------------------------------------------------
2392 : :
2393 : : !!****f* m_phonons/phonons_ncwrite
2394 : : !! NAME
2395 : : !! phonons_ncwrite
2396 : : !!
2397 : : !! FUNCTION
2398 : : !! Write phonon bandstructure to netcdf file.
2399 : : !!
2400 : : !! INPUTS
2401 : : !! ncid =NC file handle
2402 : : !! natom=Number of atoms
2403 : : !! nqpts=Number of q-points.
2404 : : !! qpoints=List of q-points in reduced coordinates
2405 : : !! weights(nqpts)= q-point weights
2406 : : !! phfreq=Phonon frequencies
2407 : : !! phdispl_cart=Phonon displacementent in Cartesian coordinates.
2408 : : !! phangmom= Phonon angular momentum in cartesian coordinates
2409 : : !!
2410 : : !! NOTES
2411 : : !! Input data is in a.u, whereas the netcdf files saves data in eV for frequencies
2412 : : !! and Angstrom for the displacements
2413 : : !! The angular momentum is output in units of hbar
2414 : : !!
2415 : : !! OUTPUT
2416 : : !! Only writing
2417 : : !!
2418 : : !! SOURCE
2419 : :
2420 : 102 : subroutine phonons_ncwrite(ncid, natom, nqpts, qpoints, weights, phfreq, phdispl_cart, phangmom)
2421 : :
2422 : : !Arguments ------------------------------------
2423 : : !scalars
2424 : : integer,intent(in) :: ncid,natom,nqpts
2425 : : !arrays
2426 : : real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts)
2427 : : real(dp),intent(in) :: phfreq(3*natom,nqpts),phdispl_cart(2,3*natom,3*natom,nqpts),phangmom(3,3*natom,nqpts)
2428 : :
2429 : : !Local variables-------------------------------
2430 : : integer :: nphmodes,ncerr
2431 : : ! *************************************************************************
2432 : :
2433 : 102 : nphmodes = 3*natom
2434 : :
2435 [ - + ]: 102 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
2436 : :
2437 : : ncerr = nctk_def_dims(ncid, [&
2438 [ + + ]: 408 : nctkdim_t("number_of_qpoints", nqpts), nctkdim_t('number_of_phonon_modes', nphmodes), nctkdim_t('three', 3)])
2439 [ - + ]: 102 : NCF_CHECK(ncerr)
2440 : :
2441 : : ! Define arrays
2442 : : ncerr = nctk_def_arrays(ncid, [&
2443 : : nctkarr_t('qpoints', "dp" , 'number_of_reduced_dimensions, number_of_qpoints'),&
2444 : : nctkarr_t('qweights',"dp", 'number_of_qpoints'),&
2445 : : nctkarr_t('phfreqs',"dp", 'number_of_phonon_modes, number_of_qpoints'),&
2446 : : nctkarr_t('phdispl_cart',"dp", 'complex, number_of_phonon_modes, number_of_phonon_modes, number_of_qpoints'),&
2447 [ + + ]: 612 : nctkarr_t('phangmom',"dp", 'three, number_of_phonon_modes, number_of_qpoints')])
2448 [ - + ]: 102 : NCF_CHECK(ncerr)
2449 : :
2450 : : ! Write variables.
2451 [ - + ]: 102 : NCF_CHECK(nctk_set_datamode(ncid))
2452 [ - + ]: 102 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints'), qpoints))
2453 [ - + ]: 102 : NCF_CHECK(nf90_put_var(ncid, vid('qweights'), weights))
2454 [ + + + + : 27523 : NCF_CHECK(nf90_put_var(ncid, vid('phfreqs'), phfreq*Ha_eV))
- + - - -
- ]
2455 [ + + + + : 777070 : NCF_CHECK(nf90_put_var(ncid, vid('phdispl_cart'), phdispl_cart*Bohr_Ang))
+ + + + -
+ - - - -
- - - - ]
2456 [ - + ]: 102 : NCF_CHECK(nf90_put_var(ncid, vid('phangmom'), phangmom))
2457 : :
2458 : : contains
2459 : 510 : integer function vid(vname)
2460 : : character(len=*),intent(in) :: vname
2461 [ + - + - ]: 510 : vid = nctk_idname(ncid, vname)
2462 : : end function vid
2463 : :
2464 : : end subroutine phonons_ncwrite
2465 : : !!***
2466 : :
2467 : : !----------------------------------------------------------------------
2468 : :
2469 : : !!****f* m_phonons/phonons_write_phfrq
2470 : : !! NAME
2471 : : !! phonons_write_phfrq
2472 : : !!
2473 : : !! FUNCTION
2474 : : !! Write phonon bandstructure in a text file. Fixed file name for the moment
2475 : : !!
2476 : : !! INPUTS
2477 : : !! natom=Number of atoms
2478 : : !! nqpts=Number of q-points.
2479 : : !! qpoints=List of q-points in reduced coordinates
2480 : : !! weights(nqpts)= q-point weights
2481 : : !! phfreq=Phonon frequencies
2482 : : !! phdispl_cart=Phonon displacementent in Cartesian coordinates.
2483 : : !!
2484 : : !! NOTES
2485 : : !! Input data is in a.u, output too
2486 : : !!
2487 : : !! OUTPUT
2488 : : !! Only writing
2489 : : !!
2490 : : !! SOURCE
2491 : :
2492 : 65 : subroutine phonons_write_phfrq(path,natom,nqpts,qpoints,weights,phfreq,phdispl_cart,phangmom)
2493 : :
2494 : : !Arguments ------------------------------------
2495 : : !scalars
2496 : : integer,intent(in) :: natom,nqpts
2497 : : character(len=*),intent(in) :: path
2498 : : !arrays
2499 : : real(dp),intent(in) :: qpoints(3,nqpts),weights(nqpts)
2500 : : real(dp),intent(in) :: phfreq(3*natom,nqpts)
2501 : : real(dp),intent(in) :: phdispl_cart(2,3*natom,3*natom,nqpts)
2502 : : real(dp),intent(in) :: phangmom(3,3*natom,nqpts)
2503 : :
2504 : : !Local variables-------------------------------
2505 : : !scalars
2506 : : integer :: nphmodes, iq, iunit, imod, icomp
2507 : : real(dp) :: dummy
2508 : : character(len=300) :: fmt
2509 : : character(len=500) :: msg
2510 : : ! *************************************************************************
2511 : :
2512 : 65 : nphmodes = 3*natom
2513 : :
2514 : 65 : dummy = qpoints(1,1); dummy = weights(1)
2515 : :
2516 : : ! Write phonon frequencies
2517 [ - + ]: 65 : if (open_file(strcat(path, "_PHFRQ"), msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2518 : 65 : ABI_ERROR(msg)
2519 : : end if
2520 : :
2521 : 65 : write (iunit, '(a)') '# ABINIT generated phonon band structure file. All in Ha atomic units'
2522 : 65 : write (iunit, '(a)') '# '
2523 : 65 : write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts
2524 : 65 : write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes
2525 : 65 : write (iunit, '(a)') '# '
2526 : :
2527 : 65 : write (fmt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)"
2528 [ + + ]: 1309 : do iq= 1, nqpts
2529 : 1309 : write (iunit, fmt) iq, phfreq(:,iq)
2530 : : end do
2531 : :
2532 : 65 : close(iunit)
2533 : :
2534 : : ! Does not Write phonon displacement ?
2535 : : if (.False.) then
2536 : : if (open_file(strcat(path, "_PHDISPL"), msg, unit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2537 : : ABI_ERROR(msg)
2538 : : end if
2539 : :
2540 : : write (iunit, '(a)') '# ABINIT generated phonon displacements, along points in PHFRQ file. All in Ha atomic units'
2541 : : write (iunit, '(a)') '# '
2542 : : write (iunit, '(a)') '# displacements in cartesian coordinates, Re and Im parts '
2543 : : write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts
2544 : : write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes
2545 : : write (iunit, '(a)') '# '
2546 : :
2547 : : !write (fmt,'(a,I3,a)') "( ", nphmodes, "(2E20.10,2x))"
2548 : : fmt = "(2E20.10,2x)"
2549 : :
2550 : : do iq = 1, nqpts
2551 : : write (iunit, '(a, i0)') '# iq ', iq
2552 : : do imod = 1, nphmodes
2553 : : write (iunit, '(a, i0)') '# imode ', imod
2554 : : do icomp = 1, nphmodes
2555 : : write (iunit, fmt, ADVANCE='NO') phdispl_cart(:,icomp,imod,iq)
2556 : : end do
2557 : : write (iunit, '(a)') ' '
2558 : : end do
2559 : : end do
2560 : :
2561 : : close(iunit)
2562 : : end if
2563 : :
2564 : : ! Write phonon angular momentum
2565 [ - + ]: 65 : if (open_file(strcat(path, "_PHANGMOM"), msg, unit=iunit, form="formatted", status="unknown", action="write") /= 0) then
2566 : 65 : ABI_ERROR(msg)
2567 : : end if
2568 : :
2569 : 65 : write (iunit, '(a)') '# ABINIT generated phonon angular momentum, along points in PHFRQ file. All in Ha atomic units'
2570 : 65 : write (iunit, '(a)') '# '
2571 : 65 : write (iunit, '(a)') '# angular momentum in cartesian coordinates '
2572 : 65 : write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts
2573 : 65 : write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes
2574 : 65 : write (iunit, '(a)') '# '
2575 : :
2576 : 65 : write (fmt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)"
2577 [ + + ]: 260 : do icomp = 1, 3
2578 [ + + ]: 3927 : do iq= 1, nqpts
2579 : 3927 : write (iunit, fmt) iq, phangmom(icomp,:,iq)
2580 : : end do
2581 [ + + ]: 260 : if (icomp /= 3) then
2582 : 130 : write (iunit, '(a,a)') ''
2583 : : end if
2584 : : end do
2585 : :
2586 : 65 : close(iunit)
2587 : :
2588 : 65 : end subroutine phonons_write_phfrq
2589 : : !!***
2590 : :
2591 : : !----------------------------------------------------------------------
2592 : :
2593 : : !!****f* m_phonons/phonons_write_xmgrace
2594 : : !! NAME
2595 : : !! phonons_write_xmgrace
2596 : : !!
2597 : : !! FUNCTION
2598 : : !! Write phonons bands in Xmgrace format. This routine should be called by a single processor.
2599 : : !!
2600 : : !! INPUTS
2601 : : !! filename=Filename
2602 : : !! natom=Number of atoms
2603 : : !! nqpts=Number of q-points
2604 : : !! qpts(3,nqpts)=Q-points
2605 : : !! phfreqs(3*natom,nqpts)=Phonon frequencies.
2606 : : !! [qptbounds(:,:)]=Optional argument giving the extrema of the q-path.
2607 : : !!
2608 : : !! OUTPUT
2609 : : !! Only writing
2610 : : !!
2611 : : !! SOURCE
2612 : :
2613 [ + + + - : 101 : subroutine phonons_write_xmgrace(filename, natom, nqpts, qpts, phfreqs, qptbounds)
- + ]
2614 : :
2615 : : !Arguments ------------------------------------
2616 : : !scalars
2617 : : integer,intent(in) :: natom,nqpts
2618 : : real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts)
2619 : : character(len=*),intent(in) :: filename
2620 : : !arrays
2621 : : real(dp),optional,intent(in) :: qptbounds(:,:)
2622 : :
2623 : : !Local variables-------------------------------
2624 : : !scalars
2625 : : integer :: unt,iq,nu,ii,start,nqbounds
2626 : : character(len=500) :: msg
2627 : : !arrays
2628 : : integer :: g0(3)
2629 : 101 : integer,allocatable :: bounds2qpt(:)
2630 : : ! *********************************************************************
2631 : :
2632 : 101 : nqbounds = 0
2633 [ + + ]: 101 : if (present(qptbounds)) then
2634 [ + + + - ]: 132 : if (product(shape(qptbounds)) > 0 ) then
2635 : : ! Find correspondence between qptbounds and k-points in ebands.
2636 : 44 : nqbounds = size(qptbounds, dim=2)
2637 [ - + + - : 132 : ABI_MALLOC(bounds2qpt, (nqbounds))
+ - + - -
+ - + ]
2638 [ + + ]: 211 : bounds2qpt = 1; start = 1
2639 [ + + ]: 211 : do ii=1,nqbounds
2640 [ + - ]: 2570 : do iq=start,nqpts
2641 [ - + - - : 2526 : if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then
- - - + +
+ ]
2642 : 167 : bounds2qpt(ii) = iq; start = iq + 1; exit
2643 : : end if
2644 : : end do
2645 : : end do
2646 : : end if
2647 : : end if
2648 : :
2649 [ - + ]: 101 : if (open_file(filename, msg, newunit=unt, form="formatted", action="write") /= 0) then
2650 : 0 : ABI_ERROR(msg)
2651 : : end if
2652 : :
2653 : 101 : write(unt,'(a)') "# Grace project file"
2654 : 101 : write(unt,'(a)') "# Generated by Abinit"
2655 : 101 : write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts
2656 : 101 : write(unt,'(a)') "# Frequencies are in meV"
2657 : 101 : write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)"
2658 [ + + ]: 3238 : do iq=1,nqpts
2659 : 3238 : write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq)))
2660 : : end do
2661 : :
2662 : 101 : write(unt,'(a)') "@page size 792, 612"
2663 : 101 : write(unt,'(a)') "@page scroll 5%"
2664 : 101 : write(unt,'(a)') "@page inout 5%"
2665 : 101 : write(unt,'(a)') "@link page off"
2666 : 101 : write(unt,'(a)') "@with g0"
2667 : 101 : write(unt,'(a)') "@world xmin 0.00"
2668 : 101 : write(unt,'(a,i0)') '@world xmax ',nqpts
2669 [ + + + + : 27358 : write(unt,'(a,e16.8)') '@world ymin ',minval(phfreqs * Ha_meV)
+ + + - -
+ - - ]
2670 [ + + + + : 27358 : write(unt,'(a,e16.8)') '@world ymax ',maxval(phfreqs * Ha_meV)
+ + + - -
+ - - ]
2671 : 101 : write(unt,'(a)') '@default linewidth 1.5'
2672 : 101 : write(unt,'(a)') '@xaxis tick on'
2673 : 101 : write(unt,'(a)') '@xaxis tick major 1'
2674 : 101 : write(unt,'(a)') '@xaxis tick major color 1'
2675 : 101 : write(unt,'(a)') '@xaxis tick major linestyle 3'
2676 : 101 : write(unt,'(a)') '@xaxis tick major grid on'
2677 : 101 : write(unt,'(a)') '@xaxis tick spec type both'
2678 : 101 : write(unt,'(a)') '@xaxis tick major 0, 0'
2679 [ + + ]: 101 : if (nqbounds /= 0) then
2680 : 44 : write(unt,'(a,i0)') '@xaxis tick spec ',nqbounds
2681 [ + + ]: 211 : do iq=1,nqbounds
2682 : : !write(unt,'(a,i0,a,a)') '@xaxis ticklabel ',iq-1,',', "foo"
2683 : 211 : write(unt,'(a,i0,a,i0)') '@xaxis tick major ',iq-1,' , ',bounds2qpt(iq) - 1
2684 : : end do
2685 : : end if
2686 : 101 : write(unt,'(a)') '@xaxis ticklabel char size 1.500000'
2687 : 101 : write(unt,'(a)') '@yaxis tick major 10'
2688 : 101 : write(unt,'(a)') '@yaxis label "Phonon Energy [meV]"'
2689 : 101 : write(unt,'(a)') '@yaxis label char size 1.500000'
2690 : 101 : write(unt,'(a)') '@yaxis ticklabel char size 1.500000'
2691 [ + + ]: 920 : do nu=1,3*natom
2692 : 920 : write(unt,'(a,i0,a)') '@ s',nu-1,' line color 1'
2693 : : end do
2694 [ + + ]: 920 : do nu=1,3*natom
2695 : 819 : write(unt,'(a,i0)') '@target G0.S',nu-1
2696 : 819 : write(unt,'(a)') '@type xy'
2697 [ + + ]: 24939 : do iq=1,nqpts
2698 : 24939 : write(unt,'(i0,1x,e16.8)') iq-1, phfreqs(nu, iq) * Ha_meV
2699 : : end do
2700 : 920 : write(unt,'(a)') '&'
2701 : : end do
2702 : :
2703 : 101 : close(unt)
2704 : :
2705 [ + + ]: 101 : ABI_SFREE(bounds2qpt)
2706 : :
2707 : 101 : end subroutine phonons_write_xmgrace
2708 : : !!***
2709 : :
2710 : : !----------------------------------------------------------------------
2711 : :
2712 : : !!****f* m_phonons/phonons_write_gnuplot
2713 : : !! NAME
2714 : : !! phonons_write_gnuplot
2715 : : !!
2716 : : !! FUNCTION
2717 : : !! Write phonons bands in gnuplot format. This routine should be called by a single processor.
2718 : : !!
2719 : : !! INPUTS
2720 : : !! prefix=prefix for files (.data, .gnuplot)
2721 : : !! natom=Number of atoms
2722 : : !! nqpts=Number of q-points
2723 : : !! qpts(3,nqpts)=Q-points
2724 : : !! phfreqs(3*natom,nqpts)=Phonon frequencies.
2725 : : !! [qptbounds(:,:)]=Optional argument giving the extrema of the q-path.
2726 : : !!
2727 : : !! OUTPUT
2728 : : !! Only writing
2729 : : !!
2730 : : !! SOURCE
2731 : :
2732 [ + - + - : 1 : subroutine phonons_write_gnuplot(prefix, natom, nqpts, qpts, phfreqs, qptbounds)
- + ]
2733 : :
2734 : : !Arguments ------------------------------------
2735 : : !scalars
2736 : : integer,intent(in) :: natom,nqpts
2737 : : real(dp),intent(in) :: qpts(3,nqpts),phfreqs(3*natom,nqpts)
2738 : : character(len=*),intent(in) :: prefix
2739 : : !arrays
2740 : : real(dp),optional,intent(in) :: qptbounds(:,:)
2741 : :
2742 : : !Local variables-------------------------------
2743 : : !scalars
2744 : : integer :: unt,iq,ii,start,nqbounds,gpl_unt
2745 : : character(len=500) :: msg,fmt
2746 : : character(len=fnlen) :: datafile,basefile
2747 : : !arrays
2748 : : integer :: g0(3)
2749 : 1 : integer,allocatable :: bounds2qpt(:)
2750 : : ! *********************************************************************
2751 : :
2752 : 1 : nqbounds = 0
2753 [ + - ]: 1 : if (present(qptbounds)) then
2754 [ + + + - ]: 3 : if (product(shape(qptbounds)) > 0 ) then
2755 : : ! Find correspondence between qptbounds and k-points in ebands.
2756 : 1 : nqbounds = size(qptbounds, dim=2)
2757 [ - + + - : 3 : ABI_MALLOC(bounds2qpt, (nqbounds))
+ - + - -
+ - + ]
2758 [ + + ]: 4 : bounds2qpt = 1; start = 1
2759 [ + + ]: 4 : do ii=1,nqbounds
2760 [ + - ]: 42 : do iq=start,nqpts
2761 [ - + - - : 41 : if (isamek(qpts(:, iq), qptbounds(:, ii), g0)) then
- - - + +
+ ]
2762 : 3 : bounds2qpt(ii) = iq; start = iq + 1; exit
2763 : : end if
2764 : : end do
2765 : : end do
2766 : : end if
2767 : : end if
2768 : :
2769 [ + - ]: 1 : datafile = strcat(prefix, "_PHBANDS.data")
2770 [ - + ]: 1 : if (open_file(datafile, msg, newunit=unt, form="formatted", action="write") /= 0) then
2771 : 0 : ABI_ERROR(msg)
2772 : : end if
2773 [ - + ]: 1 : if (open_file(strcat(prefix, "_PHBANDS.gnuplot"), msg, newunit=gpl_unt, form="formatted", action="write") /= 0) then
2774 : 1 : ABI_ERROR(msg)
2775 : : end if
2776 [ + - ]: 1 : basefile = basename(datafile)
2777 : :
2778 : 1 : write(unt,'(a)') "# Phonon band structure data file"
2779 : 1 : write(unt,'(a)') "# Generated by Abinit"
2780 : 1 : write(unt,'(2(a,i0))') "# natom: ",natom,", nqpt: ",nqpts
2781 : 1 : write(unt,'(a)') "# Frequencies are in meV"
2782 : 1 : write(unt,'(a)')"# List of q-points and their index (C notation i.e. count from 0)"
2783 [ + + ]: 42 : do iq=1,nqpts
2784 : 42 : write(unt, "(a)")sjoin("#", itoa(iq-1), ktoa(qpts(:,iq)))
2785 : : end do
2786 : :
2787 [ + - ]: 1 : fmt = sjoin("(i0,1x,", itoa(3*natom), "(es16.8,1x))")
2788 : 1 : write(unt,'(a)')"# [kpt-index, mode_1, mode_2 ...]"
2789 [ + + ]: 42 : do iq=1,nqpts
2790 [ + + ]: 165 : write(unt, fmt) iq-1, phfreqs(:, iq) * Ha_meV
2791 : : end do
2792 : :
2793 : : ! gnuplot script file
2794 : 1 : write(gpl_unt,'(a)') '# File to plot electron bandstructure with gnuplot'
2795 : : !write(gpl_unt,'(a)') "#set terminal postscript eps enhanced color font 'Times-Roman,26' lw 2"
2796 : 1 : write(gpl_unt,'(a)') '#use the next lines to make a nice figure for a paper'
2797 : 1 : write(gpl_unt,'(a)') '#set term postscript enhanced eps color lw 0.5 dl 0.5'
2798 : 1 : write(gpl_unt,'(a)') '#set pointsize 0.275'
2799 : 1 : write(gpl_unt,'(a)') 'set palette defined ( 0 "blue", 3 "green", 6 "yellow", 10 "red" )'
2800 : 1 : write(gpl_unt,'(a)') 'unset key'
2801 : 1 : write(gpl_unt,'(a)') '# can make pointsize smaller (~0.5). Too small and nothing is printed'
2802 : 1 : write(gpl_unt,'(a)') 'set pointsize 0.8'
2803 : 1 : write(gpl_unt,'(a)') 'set view 0,0'
2804 : 1 : write(gpl_unt,'(a,i0,a)') 'set xrange [0:',nqpts-1,']'
2805 : : write(gpl_unt,'(2(a,es16.8),a)')&
2806 [ + + + + : 329 : 'set yrange [',minval(phfreqs * Ha_meV),':',maxval(phfreqs * Ha_meV),']'
+ + + - -
+ - - + +
+ + + + +
- - + -
- ]
2807 : 1 : write(gpl_unt,'(a)') 'set xlabel "Momentum"'
2808 : 1 : write(gpl_unt,'(a)') 'set ylabel "Energy [meV]"'
2809 : 1 : write(gpl_unt,'(a)') strcat('set title "', replace(basefile, "_", "\\_"),'"')
2810 [ - + ]: 1 : if (nqbounds == 0) then
2811 : 0 : write(gpl_unt,'(a)') 'set grid xtics'
2812 : : else
2813 : 1 : write(gpl_unt,"(a)")"# Add vertical lines in correspondence of high-symmetry points."
2814 : 1 : write(gpl_unt,'(a)') 'unset xtics'
2815 [ + + ]: 4 : do ii=1,nqbounds
2816 : : write(gpl_unt,"(a,2(i0,a))") &
2817 : 4 : "set arrow from ",bounds2qpt(ii)-1,",graph(0,0) to ",bounds2qpt(ii)-1,",graph(1,1) nohead"
2818 : : !write(gpl_unt,"(a)")sjoin("set xtics add ('kname'", itoa(bounds2kpt(ii)-1), ")")
2819 : : end do
2820 : : end if
2821 : 1 : write(gpl_unt,"(a)")sjoin("nbranch =", itoa(3*natom))
2822 : 1 : write(gpl_unt,"(a)")strcat('plot for [i=2:nbranch] "', basefile, '" u 1:i every :1 with lines linetype -1')
2823 : 1 : write(gpl_unt,"(a)")"pause -1"
2824 : :
2825 : 1 : close(unt)
2826 : 1 : close(gpl_unt)
2827 : :
2828 [ + - ]: 1 : ABI_SFREE(bounds2qpt)
2829 : :
2830 : 1 : end subroutine phonons_write_gnuplot
2831 : : !!***
2832 : :
2833 : : !!****f* m_phonons/ifc_mkphbs
2834 : : !! NAME
2835 : : !! ifc_mkphbs
2836 : : !!
2837 : : !! FUNCTION
2838 : : !! Compute the phonon band structure from the IFC and write data to file(s)
2839 : : !!
2840 : : !! INPUTS
2841 : : !! ifc<ifc_type>=Interatomic force constants
2842 : : !! cryst<crystal_t> = Info on the crystalline structure.
2843 : : !! dtset=<datasets_type>: input: all input variables initialized from the input file.
2844 : : !! prefix=Prefix for output files.
2845 : : !! comm=MPI communicator
2846 : : !!
2847 : : !! OUTPUT
2848 : : !! Only writing.
2849 : : !!
2850 : : !! SOURCE
2851 : :
2852 : 98 : subroutine ifc_mkphbs(ifc, cryst, dtset, prefix, comm)
2853 : :
2854 : : !Arguments -------------------------------
2855 : : !scalars
2856 : : class(ifc_type),intent(in) :: ifc
2857 : : type(crystal_t),intent(in) :: cryst
2858 : : type(dataset_type),intent(in) :: dtset
2859 : : character(len=*),intent(in) :: prefix
2860 : : integer,intent(in) :: comm
2861 : :
2862 : : !Local variables -------------------------
2863 : : !scalars
2864 : : integer,parameter :: master = 0
2865 : : integer :: iqpt, nqpts, natom, ncid, nprocs, my_rank, ierr, ndirs, ncerr
2866 : 98 : type(kpath_t) :: qpath
2867 : : !arrays
2868 : 98 : real(dp),allocatable :: qph2l(:,:), qnrml2(:), eigvec(:,:,:,:,:),phfrqs(:,:),phdispl_cart(:,:,:,:),phangmom(:,:,:),weights(:)
2869 : : ! *********************************************************************
2870 : :
2871 [ - + ]: 98 : if (dtset%prtphbands == 0) return
2872 : :
2873 [ + + - + ]: 98 : if (dtset%ph_nqpath <= 0 .or. dtset%ph_ndivsm <= 0) then
2874 : 61 : ABI_COMMENT("ph_nqpath <= 0 or ph_ndivsm <= 0. Phonon bands won't be produced. Returning")
2875 : 61 : return
2876 : : end if
2877 : :
2878 : 37 : call wrtout(std_out, " Writing phonon bands, use prtphbands 0 to disable this part")
2879 : :
2880 : 37 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2881 : :
2882 : 37 : natom = cryst%natom
2883 : 37 : call qpath%init(dtset%ph_qpath(:,1:dtset%ph_nqpath), cryst%gprimd, dtset%ph_ndivsm)
2884 : 37 : nqpts = qpath%npts
2885 : :
2886 [ + - - + : 12948 : ABI_CALLOC(phfrqs, (3*natom,nqpts))
+ - - + +
- + - + -
- + - + +
+ + + ]
2887 [ - + + - : 201931 : ABI_CALLOC(phdispl_cart, (2,3*natom,3*natom,nqpts))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + ]
2888 [ - + + - : 45546 : ABI_CALLOC(phangmom, (3,3*natom,nqpts))
- + + - +
- + - - +
- + + + +
+ + + ]
2889 [ - + + - : 4357 : ABI_CALLOC(eigvec, (2,3,natom,3,natom))
+ - - + +
- + - + -
- + - + +
+ + + + +
+ + + + ]
2890 : :
2891 [ + + ]: 1971 : do iqpt=1,nqpts
2892 [ - + ]: 1934 : if (mod(iqpt, nprocs) /= my_rank) cycle ! MPI-parallelism
2893 : : ! Get phonon frequencies and displacements in cartesian coordinates for this q-point
2894 : 1934 : call ifc%fourq(cryst, qpath%points(:,iqpt), phfrqs(:,iqpt), phdispl_cart(:,:,:,iqpt), out_eigvec=eigvec)
2895 : 1971 : call phangmom_from_eigvec(natom, eigvec, phangmom(:,:,iqpt))
2896 : : end do
2897 : :
2898 : 37 : call xmpi_sum_master(phfrqs, master, comm, ierr)
2899 : 37 : call xmpi_sum_master(phdispl_cart, master, comm, ierr)
2900 : 37 : call xmpi_sum_master(phangmom, master, comm, ierr)
2901 : :
2902 [ + - ]: 37 : if (my_rank == master) then
2903 [ - + + - : 111 : ABI_MALLOC(weights, (nqpts))
+ - + - -
+ - + ]
2904 [ + + ]: 1971 : weights = one
2905 : :
2906 : : ! Compute directions for non-analytical behaviour.
2907 : : ! TODO: The same approach should be used in anaddb at the level of the parser.
2908 [ - + + - : 111 : ABI_MALLOC(qph2l, (3, 2*dtset%ph_nqpath))
+ - - + -
+ ]
2909 [ - + + - : 111 : ABI_MALLOC(qnrml2, (2*dtset%ph_nqpath))
+ - - + -
+ ]
2910 : :
2911 : 37 : ndirs = 0
2912 [ + + + + : 531 : if (any(ifc%zeff /= zero)) then
+ + + + +
+ ]
2913 [ + + ]: 73 : do iqpt=1,dtset%ph_nqpath
2914 [ + + + + ]: 247 : if (sum(dtset%ph_qpath(:, iqpt)**2) < tol14) then
2915 : 16 : ndirs = ndirs + 1
2916 [ + + ]: 16 : if (iqpt == 1) then
2917 [ + + ]: 60 : qph2l(:, ndirs) = dtset%ph_qpath(:, 2) - dtset%ph_qpath(:, 1)
2918 [ - + ]: 1 : else if (iqpt == dtset%ph_nqpath) then
2919 [ # # ]: 0 : qph2l(:, ndirs) = dtset%ph_qpath(:, dtset%ph_nqpath - 1) - dtset%ph_qpath(:, dtset%ph_nqpath)
2920 : : else
2921 [ + + ]: 4 : qph2l(:, ndirs) = dtset%ph_qpath(:, iqpt - 1) - dtset%ph_qpath(:, iqpt)
2922 : 1 : ndirs = ndirs + 1
2923 [ + + ]: 4 : qph2l(:, ndirs) = dtset%ph_qpath(:, iqpt + 1) - dtset%ph_qpath(:, iqpt)
2924 : : end if
2925 : : end if
2926 : : end do
2927 : :
2928 : : ! Convert to Cartesian coordinates.
2929 [ + + ]: 32 : do iqpt=1,ndirs
2930 [ + + + + : 338 : qph2l(:, iqpt) = matmul(cryst%gprimd, qph2l(:, iqpt))
+ + + + +
- ]
2931 : : end do
2932 [ + + ]: 131 : qnrml2 = zero
2933 : : end if
2934 : :
2935 : : ! TODO: A similar piece of code is used in anaddb (mkpbs + ifc_calcnwrite_nana_terms).
2936 : : ! Should centralize everything in a single routine
2937 [ - + ]: 37 : NCF_CHECK_MSG(nctk_open_create(ncid, strcat(prefix, "_PHBST.nc"), xmpi_comm_self), "Creating PHBST")
2938 [ - + ]: 37 : NCF_CHECK(cryst%ncwrite(ncid))
2939 : 37 : call phonons_ncwrite(ncid, natom, nqpts, qpath%points, weights, phfrqs, phdispl_cart, phangmom)
2940 : : ! This flag tells AbiPy that all the non-analytic directions have been computed.
2941 [ + + - + : 74 : NCF_CHECK(nctk_defnwrite_ivars(ncid, ["has_abipy_non_anal_ph"], [1]))
- - ]
2942 [ + + ]: 74 : ncerr = nctk_def_arrays(ncid, [nctkarr_t("atomic_mass_units", "dp", "number_of_atom_species")], defmode=.True.)
2943 [ - + ]: 37 : NCF_CHECK(ncerr)
2944 [ - + ]: 37 : NCF_CHECK(nctk_set_datamode(ncid))
2945 [ - + ]: 37 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "atomic_mass_units"), ifc%amu))
2946 [ + + ]: 37 : if (ndirs /= 0) call ifc%calcnwrite_nana_terms(cryst, ndirs, qph2l, qnrml2, ncid=ncid)
2947 [ - + ]: 37 : NCF_CHECK(nf90_close(ncid))
2948 : :
2949 [ - + ]: 37 : ABI_FREE(qph2l)
2950 [ - + ]: 37 : ABI_FREE(qnrml2)
2951 : :
2952 : 36 : select case (dtset%prtphbands)
2953 : : case (1)
2954 : 36 : call phonons_write_xmgrace(strcat(prefix, "_PHBANDS.agr"), natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds)
2955 : : case (2)
2956 : 1 : call phonons_write_gnuplot(prefix, natom, nqpts, qpath%points, phfrqs, qptbounds=qpath%bounds)
2957 : : case (3)
2958 : 0 : call phonons_write_phfrq(prefix, natom, nqpts, qpath%points, weights, phfrqs, phdispl_cart, phangmom)
2959 : : case default
2960 [ + + - - ]: 37 : ABI_WARNING(sjoin("Unsupported value for prtphbands:", itoa(dtset%prtphbands)))
2961 : : end select
2962 : :
2963 [ - + ]: 37 : ABI_FREE(weights)
2964 : : end if ! master
2965 : :
2966 [ - + ]: 37 : ABI_FREE(phfrqs)
2967 [ - + ]: 37 : ABI_FREE(phdispl_cart)
2968 [ - + ]: 37 : ABI_FREE(phangmom)
2969 [ - + ]: 37 : ABI_FREE(eigvec)
2970 : :
2971 : 37 : call qpath%free()
2972 : :
2973 [ - + - + : 98 : end subroutine ifc_mkphbs
- + - + -
+ - + - +
- + ]
2974 : : !!***
2975 : :
2976 : : !!****f* m_phonons/dfpt_symph
2977 : : !! NAME
2978 : : !! dfpt_symph
2979 : : !!
2980 : : !! FUNCTION
2981 : : !! Determine the symmetry character of the different phonon modes.
2982 : : !!
2983 : : !! INPUTS
2984 : : !! acell(3)=length scales of primitive translations (bohr)
2985 : : !! eigvec(2*3*natom*3*natom)=eigenvectors of the dynamical matrix
2986 : : !! indsym(4,nsym,natom)=indirect indexing array : for each
2987 : : !! isym,iatom, fourth element is label of atom into which iatom is sent by
2988 : : !! INVERSE of symmetry operation isym; first three elements are the primitive
2989 : : !! translations which must be subtracted after the transformation to get back
2990 : : !! to the original unit cell.
2991 : : !! iout=unit number to which output is written
2992 : : !! natom=number of atoms in unit cell
2993 : : !! nsym=number of space group symmetries
2994 : : !! phfrq(3*natom)=phonon frequencies (Hartree units)
2995 : : !! rprim(3,3)=dimensionless primitive translations in real space
2996 : : !! symrel(3,3,nsym)=matrices of the group symmetries (real space)
2997 : : !!
2998 : : !! OUTPUT
2999 : : !!
3000 : : !! SOURCE
3001 : :
3002 : 72 : subroutine dfpt_symph(iout, acell, eigvec, indsym, natom, nsym, phfrq, rprim, symrel)
3003 : :
3004 : : !Arguments ------------------------------------
3005 : : !scalars
3006 : : integer,intent(in) :: iout,natom,nsym
3007 : : !arrays
3008 : : integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym)
3009 : : real(dp),intent(in) :: acell(3),eigvec(2*3*natom*3*natom),phfrq(3*natom)
3010 : : real(dp),intent(in) :: rprim(3,3)
3011 : :
3012 : : !Local variables -------------------------
3013 : : !scalars
3014 : : integer :: iad1,iad2,iad3,iatom,idir,ii1,ii2,ii3,imode,isym,itol,jad,jatom,jj, jmode,kk,ntol
3015 : : character(len=500) :: msg
3016 : : !arrays
3017 : : integer :: units(2)
3018 : 72 : integer,allocatable :: degeneracy(:),integer_characters(:),symind(:,:)
3019 : : real(dp) :: gprimd(3,3),rprimd(3,3)
3020 : 72 : real(dp),allocatable :: eigvtr(:),redvec(:),redvtr(:),symph(:,:)
3021 : : !******************************************************************
3022 : :
3023 [ + + ]: 216 : units = [std_out, iout]
3024 : :
3025 : : ! Compute dimensional primitive translations rprimd and its inverse gprimd
3026 : 72 : call mkrdim(acell,rprim,rprimd)
3027 : 72 : call matr3inv(rprimd,gprimd)
3028 : :
3029 : : ! Build the symmetry index (inverse of indsym(4,:,:))
3030 [ - + + - : 288 : ABI_MALLOC(symind, (nsym,natom))
- + + - +
- + - - +
- + ]
3031 [ + + ]: 2275 : do isym=1,nsym
3032 [ + + ]: 7984 : do iatom=1,natom
3033 : 7912 : symind(isym,indsym(4,isym,iatom))=iatom
3034 : : end do
3035 : : end do
3036 : :
3037 [ - + + - : 360 : ABI_MALLOC(symph,(nsym,3*natom))
+ - - + +
- + - + -
- + - + ]
3038 [ + - - + : 216 : ABI_MALLOC(redvec,(2*3*natom))
+ - + - +
- - + -
+ ]
3039 [ - + + - : 144 : ABI_MALLOC(redvtr,(2*3*natom))
- + - + ]
3040 [ - + + - : 144 : ABI_MALLOC(eigvtr,(2*3*natom))
- + - + ]
3041 : :
3042 : : ! Loop over the vibration modes
3043 [ + + ]: 756 : do imode=1,3*natom
3044 : :
3045 : : ! Compute eigvec for this mode in reduced coordinates redvec
3046 [ + + ]: 4320 : do iatom=1,natom
3047 : 3636 : iad1=3*(iatom-1)+1
3048 : 3636 : ii1=2*3*natom*(imode-1)+2*(iad1-1)+1
3049 : 3636 : iad2=3*(iatom-1)+2
3050 : 3636 : ii2=2*3*natom*(imode-1)+2*(iad2-1)+1
3051 : 3636 : iad3=3*(iatom-1)+3
3052 : 3636 : ii3=2*3*natom*(imode-1)+2*(iad3-1)+1
3053 [ + + ]: 15228 : do idir=1,3
3054 : 10908 : jad=3*(iatom-1)+idir
3055 : 10908 : jj=2*(jad-1)+1
3056 : : redvec(jj)=gprimd(1,idir)*eigvec(ii1)+&
3057 : : gprimd(2,idir)*eigvec(ii2)+&
3058 : 10908 : gprimd(3,idir)*eigvec(ii3)
3059 : : redvec(jj+1)=gprimd(1,idir)*eigvec(ii1+1)+&
3060 : : gprimd(2,idir)*eigvec(ii2+1)+&
3061 : 14544 : gprimd(3,idir)*eigvec(ii3+1)
3062 : : end do !idir
3063 : : end do !iatom
3064 : :
3065 : : ! Apply each transformation to redvec and store at the correct location in redvtr (iatom -> jatom)
3066 [ + + ]: 17883 : do isym=1,nsym
3067 [ + + ]: 90960 : do iatom=1,natom
3068 : 73833 : jatom=symind(isym,iatom)
3069 : 73833 : iad1=3*(iatom-1)+1
3070 : 73833 : ii1=2*(iad1-1)+1
3071 : 73833 : iad2=3*(iatom-1)+2
3072 : 73833 : ii2=2*(iad2-1)+1
3073 : 73833 : iad3=3*(iatom-1)+3
3074 : 73833 : ii3=2*(iad3-1)+1
3075 [ + + ]: 312459 : do idir=1,3
3076 : 221499 : jad=3*(jatom-1)+idir
3077 : 221499 : jj=2*(jad-1)+1
3078 : : redvtr(jj)=dble(symrel(idir,1,isym))*redvec(ii1)+&
3079 : : dble(symrel(idir,2,isym))*redvec(ii2)+&
3080 : 221499 : dble(symrel(idir,3,isym))*redvec(ii3)
3081 : : redvtr(jj+1)=dble(symrel(idir,1,isym))*redvec(ii1+1)+&
3082 : : dble(symrel(idir,2,isym))*redvec(ii2+1)+&
3083 : 295332 : dble(symrel(idir,3,isym))*redvec(ii3+1)
3084 : : end do !idir
3085 : : end do !iatom
3086 : :
3087 : : ! Compute redvtr in cartesian coordinates eigvtr
3088 [ + + ]: 90960 : do iatom=1,natom
3089 : 73833 : iad1=3*(iatom-1)+1
3090 : 73833 : ii1=2*(iad1-1)+1
3091 : 73833 : iad2=3*(iatom-1)+2
3092 : 73833 : ii2=2*(iad2-1)+1
3093 : 73833 : iad3=3*(iatom-1)+3
3094 : 73833 : ii3=2*(iad3-1)+1
3095 [ + + ]: 312459 : do idir=1,3
3096 : 221499 : jad=3*(iatom-1)+idir
3097 : 221499 : jj=2*(jad-1)+1
3098 : : eigvtr(jj)=rprimd(idir,1)*redvtr(ii1)+&
3099 : : rprimd(idir,2)*redvtr(ii2)+&
3100 : 221499 : rprimd(idir,3)*redvtr(ii3)
3101 : : eigvtr(jj+1)=rprimd(idir,1)*redvtr(ii1+1)+&
3102 : : rprimd(idir,2)*redvtr(ii2+1)+&
3103 : 295332 : rprimd(idir,3)*redvtr(ii3+1)
3104 : : end do !idir
3105 : : end do !iatom
3106 : :
3107 : : ! Compute scalar product...
3108 : 17127 : symph(isym,imode)=zero
3109 [ + + ]: 239310 : do jad=1,3*natom
3110 : 221499 : jj=2*(jad-1)+1
3111 : 221499 : kk=2*3*natom*(imode-1)+2*(jad-1)+1
3112 : 238626 : symph(isym,imode)=symph(isym,imode)+eigvtr(jj)*eigvec(kk)+eigvtr(jj+1)*eigvec(kk+1)
3113 : : end do
3114 : :
3115 : : end do !isym
3116 : : end do !imode
3117 : :
3118 : : !Treat degeneracies (different tolerances will be tried)
3119 : : !Compute the order of the degeneracy, and
3120 : : !attribute it to the lowest of the degenerate modes
3121 : : !Also attribute the characters to the lowest mode
3122 : : !When all the characters are integers, consider that the
3123 : : !mode is non-degenerate. The maximum difference in frequency
3124 : : !that is tolerated is on the order of 4cm-1 (which is large...)
3125 [ - + + - : 216 : ABI_MALLOC(degeneracy, (3*natom))
+ - + - -
+ - + ]
3126 [ - + + - : 144 : ABI_MALLOC(integer_characters, (3*natom))
- + - + ]
3127 [ + + ]: 756 : degeneracy(:)=1
3128 [ + + ]: 756 : integer_characters(:)=0
3129 [ + + ]: 121 : do itol=1,20
3130 : 119 : ntol=itol
3131 [ + + ]: 1155 : do imode=3*natom,2,-1
3132 [ + + ]: 1155 : if(integer_characters(imode)==0)then
3133 [ + + ]: 6187 : do jmode=imode-1,1,-1
3134 [ + + ]: 6187 : if(integer_characters(jmode)==0)then
3135 [ + + ]: 5199 : if(abs(phfrq(imode)-phfrq(jmode))<itol*tol6)then
3136 : 523 : degeneracy(jmode)=degeneracy(jmode)+degeneracy(imode)
3137 : 523 : degeneracy(imode)=0
3138 [ + + ]: 16241 : symph(:,jmode)=symph(:,jmode)+symph(:,imode)
3139 [ + + ]: 16241 : symph(:,imode)=zero
3140 : : end if
3141 : : end if !integer_characters(jmode)==0
3142 : : end do !jmode
3143 : : end if !integer_characters(imode)==0
3144 : : end do !imode
3145 [ + + ]: 1274 : do imode=1,3*natom
3146 [ + - - + : 22838 : if(maxval(abs( symph(:,imode)-nint(symph(:,imode)) ))<0.05_dp)then
- - + + +
+ ]
3147 : 1041 : integer_characters(imode)=1
3148 : : end if
3149 : : end do
3150 [ + + + + ]: 1276 : if(sum(integer_characters(:))==3*natom)exit
3151 : : end do !itol
3152 : :
3153 : : !write(std_out,*)' dfpt_symph : degeneracy=',degeneracy(:)
3154 : :
3155 : 72 : write(msg,'(a,a,es8.2,5a)')ch10,' Analysis of degeneracies and characters (maximum tolerance=',ntol*tol6,' a.u.)',ch10,&
3156 : 72 : ' For each vibration mode, or group of modes if degenerate,',ch10,&
3157 : 144 : ' the characters are given for each symmetry operation (see the list in the log file).'
3158 : 72 : call wrtout(units, msg)
3159 : :
3160 : : !use m_ptgroups, only : get_classes
3161 : : !call get_classes(nsym, sym, nclass, nelements, elements_idx)
3162 : :
3163 [ + + ]: 756 : do imode=1,3*natom
3164 [ + + ]: 756 : if(degeneracy(imode)/=0)then
3165 : 334 : write(msg,'(a,i4)') ' Symmetry characters of vibration mode #',imode
3166 : 334 : call wrtout(units, msg)
3167 [ + + ]: 334 : if (degeneracy(imode)>=2)then
3168 [ + + ]: 201 : if (degeneracy(imode)==2) write(msg,'(a,i4)') ' degenerate with vibration mode #',imode+1
3169 [ + + ]: 201 : if (degeneracy(imode)>=3) write(msg,'(a,i4,a,i4)') &
3170 : 137 : ' degenerate with vibration modes #',imode+1,' to ',imode+degeneracy(imode)-1
3171 : 201 : call wrtout(units, msg)
3172 : : end if
3173 [ + + ]: 880 : do jj=1,(nsym-1)/16+1
3174 : 546 : write(msg,'(16f5.1)') (symph(isym,imode),isym=(jj-1)*16+1,min(nsym,jj*16))
3175 : 880 : call wrtout(units, msg)
3176 : : end do
3177 : : end if
3178 : : end do !imode
3179 : :
3180 [ - + ]: 72 : ABI_FREE(degeneracy)
3181 [ - + ]: 72 : ABI_FREE(integer_characters)
3182 [ - + ]: 72 : ABI_FREE(eigvtr)
3183 [ - + ]: 72 : ABI_FREE(redvtr)
3184 [ - + ]: 72 : ABI_FREE(redvec)
3185 [ - + ]: 72 : ABI_FREE(symph)
3186 [ - + ]: 72 : ABI_FREE(symind)
3187 : :
3188 : 72 : end subroutine dfpt_symph
3189 : : !!***
3190 : :
3191 : : !!****f* m_phonons/freeze_displ_allmodes
3192 : : !!
3193 : : !! NAME
3194 : : !! freeze_displ_allmodes
3195 : : !!
3196 : : !! FUNCTION
3197 : : !! From a given set of phonon modes, generate and output supercells and
3198 : : !! displaced configurations of atoms.
3199 : : !! Typically useful to follow soft modes and see distortions of crystal structures
3200 : : !!
3201 : : !! INPUTS
3202 : : !! amu(ntypat) = mass of the atoms (atomic mass unit)
3203 : : !! displ(2,3*natom,3*natom) = phonon mode displacements (complex)
3204 : : !! freeze_displ = amplitude of the displacement to freeze into the supercell
3205 : : !! natom = number of atoms in the unit cell
3206 : : !! ntypat = number of atom types
3207 : : !! phfrq(3*natom) = phonon frequencies
3208 : : !! qphnrm = norm of phonon q vector (should be 1 or 0)
3209 : : !! qphon = phonon wavevector
3210 : : !! rprimd(3,3) = dimensionfull primitive translations in real space
3211 : : !! typat(natom) = integer label of each type of atom (1,2,...)
3212 : : !! xcart(3,natom) = cartesian coords of atoms in unit cell (bohr)
3213 : : !!
3214 : : !! OUTPUT
3215 : : !! for the moment only prints to file, but could also return pointer to supercell object, with
3216 : : !! rprimd and atomic positions, for further use
3217 : : !!
3218 : : !! NOTES
3219 : : !! freeze_displ could be determined automatically from a temperature and the phonon frequency,
3220 : : !! as the average displacement of the mode with a Bose distribution.
3221 : : !!
3222 : : !! SOURCE
3223 : : !!
3224 : :
3225 : 12 : subroutine freeze_displ_allmodes(displ, freeze_displ, natom, outfile_radix, phfreq, &
3226 [ - + ]: 6 : qphon, rprimd, typat, xcart, znucl)
3227 : :
3228 : : !Arguments ------------------------------------
3229 : : !scalars
3230 : : integer,intent(in) :: natom
3231 : : character(len=*),intent(in) :: outfile_radix
3232 : : real(dp), intent(in) :: freeze_displ
3233 : : !arrays
3234 : : integer,intent(in) :: typat(natom)
3235 : : real(dp),intent(in) :: displ(2,3*natom,3*natom), rprimd(3,3), phfreq(3*natom), qphon(3), xcart(3,natom)
3236 : : real(dp),intent(in) :: znucl(:)
3237 : :
3238 : : !Local variables-------------------------------
3239 : : integer :: jmode
3240 : 6 : type(supercell_type) :: scell
3241 : : ! *************************************************************************
3242 : :
3243 : : !determine supercell needed to freeze phonon
3244 : 6 : call scell%init_for_qpt(natom, qphon, rprimd, typat, xcart, znucl)
3245 : :
3246 [ + + ]: 42 : do jmode = 1, 3*natom
3247 : : ! reset positions
3248 [ + - + - : 504 : scell%xcart = scell%xcart_ref
- + - - -
- - - + +
+ + ]
3249 : :
3250 : : ! displace atoms according to phonon jmode
3251 : 36 : call scell%freeze_displ(displ(:,:,jmode), freeze_displ)
3252 : :
3253 : : ! print out everything for this wavevector and mode
3254 : 42 : call scell%print_for_qpt(phfreq(jmode), jmode, outfile_radix)
3255 : : end do
3256 : :
3257 : 6 : call scell%free()
3258 : :
3259 [ - + - + : 6 : end subroutine freeze_displ_allmodes
- + - + -
+ - + -
+ ]
3260 : : !!***
3261 : :
3262 : : !----------------------------------------------------------------------
3263 : :
3264 : : !!****f* m_phonons/pheigvec_rotate
3265 : : !! NAME
3266 : : !! pheigvec_rotate
3267 : : !!
3268 : : !! FUNCTION
3269 : : !! Return phonon eigenvectors for q in the BZ from the symmetrical image in the IBZ.
3270 : : !! Note that the isym index is supposed to have been computed in kpts_map with the symrec convention.
3271 : : !!
3272 : : !! INPUTS
3273 : : !! cryst: crystal structure
3274 : : !! qq_ibz: q-point in the IBZ
3275 : : !! isym
3276 : : !! itimrev
3277 : : !! eigvec_ibz: Input phonon eigenvectors at qq_ibz.
3278 : : !!
3279 : : !! OUTPUT
3280 : : !! eigvec_bz: phonon eigenvectors at q_bz.
3281 : : !! displ_cart_qbz: phonon displacement at q_bz in Cartesian coordinates.
3282 : : !! [displ_red_qbz]: phonon displacement at q_bz in reduced coordinates.
3283 : :
3284 : 4777 : subroutine pheigvec_rotate(cryst, qq_ibz, isym, itimrev, eigvec_ibz, eigvec_qbz, displ_cart_qbz, &
3285 : : displ_red_qbz) ! Optional
3286 : :
3287 : : !Arguments ------------------------------------
3288 : : !scalars
3289 : : type(crystal_t),intent(in) :: cryst
3290 : : integer,intent(in) :: isym, itimrev
3291 : : real(dp),intent(in) :: qq_ibz(3), eigvec_ibz(2,3*cryst%natom,3*cryst%natom)
3292 : : real(dp),intent(out) :: eigvec_qbz(2,3*cryst%natom,3*cryst%natom)
3293 : : real(dp),intent(out) :: displ_cart_qbz(2,3*cryst%natom,3*cryst%natom)
3294 : : real(dp),optional,intent(out) :: displ_red_qbz(2,3*cryst%natom,3*cryst%natom)
3295 : :
3296 : : !Local variables-------------------------------
3297 : : !scalars
3298 : : integer :: natom, natom3, idir, iat, jdir, iat_sym !, isym_inv
3299 : : real(dp) :: arg
3300 : : !arrays
3301 : : integer :: r0(3)
3302 : 9554 : real(dp) :: gamma_matrix(2,3,cryst%natom,3,cryst%natom)
3303 : : real(dp) :: symat(3,3), phase(2) !, dum(0, 0), gamma2(2,3,cryst%natom,3,cryst%natom)
3304 : : !************************************************************************
3305 : :
3306 : 4777 : natom = cryst%natom; natom3 = cryst%natom * 3
3307 [ + + + + ]: 62101 : symat = cryst%symrel_cart(:,:,isym)
3308 : :
3309 : : ! Build Gamma matrix in Cartesian coordinates.
3310 : : ! e(S q_ibz) = Gamma({S, v}) e(q_ibz)
3311 [ + + + + : 542819 : gamma_matrix = zero
+ + + + +
+ ]
3312 [ + + ]: 13550 : do iat=1,natom
3313 : : !do iat_sym=1,natom
3314 : : ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
3315 : : ! indsym(4, isym,iat) gives iat_sym in the original unit cell.
3316 : : ! indsym(1:3,isym,iat) gives the lattice vector $R_0$.
3317 : 8773 : iat_sym = cryst%indsym(4, isym, iat)
3318 [ + + ]: 35092 : r0 = cryst%indsym(1:3, isym, iat)
3319 [ + + ]: 35092 : arg = two_pi * dot_product(qq_ibz, real(r0))
3320 [ + + ]: 26319 : phase(:) = [cos(arg), sin(arg)] !; write(std_out, *)" ro: ", r0, "phase: " ,phase, "qq_ibz: ", qq_ibz
3321 [ + + ]: 39869 : do jdir=1,3
3322 [ + + ]: 114049 : do idir=1,3
3323 [ + + ]: 263190 : gamma_matrix(:, idir, iat, jdir, iat_sym) = symat(idir, jdir) * phase(:)
3324 : : end do
3325 : : end do
3326 : : end do
3327 : :
3328 : : !write(std_out, "(2a)")" Gamma_matrix for qq_bz:", trim(ktoa(qq_bz))
3329 : : !call print_arr([std_out], reshape(cmplx(gamma_matrix(1,:,:,:,:), gamma_matrix(2,:,:,:,:)), [natom3, natom3]))
3330 : : !gamma2 = gamma_matrix
3331 : : !call cg_zgemm("C", "N", natom3, natom3, natom3, gamma_matrix, gamma2, eigvec_qbz)
3332 : : !write(std_out, "(a)")" gamma^H gamma:"
3333 : : !call print_arr([std_out], reshape(cmplx(eigvec_qbz(1,:,:), eigvec_qbz(2,:,:)), [natom3, natom3]))
3334 : : !call cg_check_unitary(natom3, gamm_matrix)
3335 : :
3336 : 4777 : call cg_zgemm("N", "N", natom3, natom3, natom3, gamma_matrix, eigvec_ibz, eigvec_qbz)
3337 [ + + + + : 15781 : if (itimrev == 1) eigvec_qbz(2,:,:) = -eigvec_qbz(2,:,:)
+ + ]
3338 : :
3339 : : ! Fix the phase of the eigenvectors
3340 : : !call fxphas_seq(eigvec_qbz, dum, 0, 0, 1, 3*natom*3*natom, 0, 3*natom, 3*natom, 0)
3341 : : ! Normalise the eigenvectors
3342 : : !call pheigvec_normalize(natom, eigvec_qbz)
3343 : :
3344 : : ! Compute phonon displacements in Cartesian coordinates
3345 : 4777 : call phdispl_from_eigvec(cryst%natom, cryst%ntypat, cryst%typat, cryst%amu, eigvec_qbz, displ_cart_qbz)
3346 : :
3347 : : ! Compute phonon displacements in reduced coordinates.
3348 [ + + ]: 4777 : if (present(displ_red_qbz)) call phdispl_cart2red(cryst%natom, cryst%gprimd, displ_cart_qbz, displ_red_qbz)
3349 : :
3350 : 4777 : end subroutine pheigvec_rotate
3351 : : !!***
3352 : :
3353 : : !----------------------------------------------------------------------
3354 : :
3355 : : !!****f* m_phonons/phstore_init
3356 : : !! NAME
3357 : : !! phstore_init
3358 : : !!
3359 : : !! FUNCTION
3360 : : !! Create new object with phonon quantities in the IBZ.
3361 : : !!
3362 : : !! INPUTS
3363 : : !! cryst: Crystal structure
3364 : : !! ifc: Interatomic force constants.
3365 : : !! nqibz: Number of q-points in the IBZ.
3366 : : !! qibz: q-points in the IBZ
3367 : : !! use_ifc_fourq: True to replace symmetrization with call to ifc_fourq (debugging option)
3368 : : !! comm: MPI communicator in which phonon arrays in the IBZ will be MPI distributed.
3369 : :
3370 [ + - ]: 58 : subroutine phstore_init(new, cryst, ifc, nqibz, qibz, use_ifc_fourq, comm)
3371 : :
3372 : : !Arguments ------------------------------------
3373 : : class(phstore_t),intent(out) :: new
3374 : : type(crystal_t),intent(in) :: cryst
3375 : : type(ifc_type),intent(in) :: ifc
3376 : : integer,intent(in) :: nqibz, comm
3377 : : logical,intent(in) :: use_ifc_fourq
3378 : : real(dp),target,intent(in) :: qibz(3, nqibz)
3379 : :
3380 : : !Local variables ------------------------------
3381 : : !scalars
3382 : : integer :: natom3, my_q1, my_q2, iq_ibz
3383 : : character(len=500) :: msg
3384 : : ! *************************************************************************
3385 : :
3386 : 58 : new%qibz => qibz
3387 : :
3388 : 58 : new%natom = cryst%natom; natom3 = cryst%natom * 3; new%natom3 = natom3
3389 : 58 : new%comm = comm; new%nprocs = xmpi_comm_size(comm); new%my_rank = xmpi_comm_rank(comm)
3390 : 58 : new%use_ifc_fourq = use_ifc_fourq
3391 : :
3392 [ - + + - : 232 : ABI_MALLOC(new%displ_cart, (2, 3, cryst%natom, natom3))
- + + - +
- + - - +
- + - + ]
3393 [ + - - + : 232 : ABI_MALLOC(new%phfrq, (3*cryst%natom))
+ - + - +
- - + - +
- + ]
3394 [ - + ]: 58 : if (new%use_ifc_fourq) return
3395 : :
3396 : : ! Split qibz in blocks inside comm
3397 [ - + + - : 174 : ABI_MALLOC(new%qibz_start, (0:new%nprocs-1))
+ - + - -
+ - + -
+ ]
3398 [ - + + - : 116 : ABI_MALLOC(new%qibz_stop, (0:new%nprocs-1))
- + - + -
+ ]
3399 : 58 : call xmpi_split_work2_i4b(nqibz, new%nprocs, new%qibz_start, new%qibz_stop)
3400 : :
3401 : 58 : my_q1 = new%qibz_start(new%my_rank)
3402 : 58 : my_q2 = new%qibz_stop(new%my_rank)
3403 : :
3404 : 58 : call wrtout(std_out, " Computing all phonon frequencies and eigenvectors in the IBZ.", pre_newlines=1)
3405 : 58 : call wrtout(std_out, sjoin(" Number of IBZ q-points stored by this rank inside pert_comm:", itoa(my_q2 - my_q1 + 1)))
3406 : : write(msg, "(a,f8.1,a)") &
3407 : 58 : " Memory required by pheigvec_qibz: ", 2 * natom3**2 * (my_q2 - my_q1 + 1) * dp * b2Mb, " [Mb] <<< MEM"
3408 : 58 : call wrtout(std_out, msg)
3409 : :
3410 [ - + + - : 232 : ABI_MALLOC(new%phfreqs_qibz, (natom3, my_q1:my_q2))
- + + - +
- + - - +
- + - + ]
3411 [ - + + - : 290 : ABI_MALLOC(new%pheigvec_qibz, (2, natom3, natom3, my_q1:my_q2))
- + + - -
+ + - + -
+ - - + -
+ - + ]
3412 : :
3413 [ + + ]: 4048 : do iq_ibz=my_q1, my_q2
3414 : : call ifc%fourq(cryst, qibz(:,iq_ibz), new%phfreqs_qibz(:, iq_ibz), new%displ_cart, &
3415 : 4048 : out_eigvec=new%pheigvec_qibz(:,:,:,iq_ibz))
3416 : : end do
3417 : :
3418 : : end subroutine phstore_init
3419 : : !!***
3420 : :
3421 : : !----------------------------------------------------------------------
3422 : :
3423 : : !!****f* m_phonons/phstore_free
3424 : : !! NAME
3425 : : !! phstore_free
3426 : : !!
3427 : : !! FUNCTION
3428 : : !! Free dynamic memory.
3429 : : !!
3430 : : !! INPUTS
3431 : :
3432 : 58 : subroutine phstore_free(self)
3433 : :
3434 : : !Arguments ------------------------------------
3435 : : class(phstore_t),intent(inout) :: self
3436 : : ! *************************************************************************
3437 : :
3438 [ + - ]: 58 : ABI_SFREE(self%qibz_start)
3439 [ + - ]: 58 : ABI_SFREE(self%qibz_stop)
3440 [ + - ]: 58 : ABI_SFREE(self%phfreqs_qibz)
3441 [ + - ]: 58 : ABI_SFREE(self%pheigvec_qibz)
3442 [ + - ]: 58 : ABI_SFREE(self%displ_cart)
3443 [ + - ]: 58 : ABI_SFREE(self%phfrq)
3444 : :
3445 : 58 : self%qibz => null()
3446 : :
3447 : 58 : end subroutine phstore_free
3448 : : !!***
3449 : :
3450 : : !----------------------------------------------------------------------
3451 : :
3452 : : !!****f* m_phonons/phstore_async_rotate
3453 : : !! NAME
3454 : : !! phstore_async_rotate
3455 : : !!
3456 : : !! FUNCTION
3457 : : !! Begin non-blocking collective MPI communication inside self%comm to obtain
3458 : : !! phonon frequencies and eigenvectors in the BZ from data in the IBZ.
3459 : : !!
3460 : : !! INPUTS
3461 : :
3462 : 4852 : subroutine phstore_async_rotate(self, cryst, ifc, iq_ibz, qpt_ibz, qpt_bz, isym_q, trev_q)
3463 : :
3464 : : !Arguments ------------------------------------
3465 : : class(phstore_t),intent(inout) :: self
3466 : : type(crystal_t), intent(in) :: cryst
3467 : : type(ifc_type),intent(in) :: ifc
3468 : : integer,intent(in) :: iq_ibz, isym_q, trev_q
3469 : : real(dp),intent(in) :: qpt_ibz(3), qpt_bz(3)
3470 : :
3471 : : !Local variables ------------------------------
3472 : : !scalars
3473 : : integer :: rank, master, ierr
3474 : : logical :: isirr_q
3475 : 4852 : real(dp) :: eigvec_qpt(2, self%natom3, self%natom3)
3476 : : ! *************************************************************************
3477 : :
3478 : : ABI_UNUSED(qpt_ibz(1))
3479 : :
3480 [ - + ]: 4852 : if (self%use_ifc_fourq) then
3481 : : ! Debugging section.
3482 : 0 : call ifc%fourq(cryst, qpt_bz, self%phfrq, self%displ_cart); return
3483 : : end if
3484 : :
3485 : : ! Find the MPI rank storing the q-point in the IBZ.
3486 [ + - ]: 4852 : do rank=0,self%nprocs-1
3487 [ + - + - ]: 4852 : if (iq_ibz >= self%qibz_start(rank) .and. iq_ibz <= self%qibz_stop(rank)) then
3488 : 4852 : master = rank; exit
3489 : : end if
3490 : : end do
3491 [ - + ]: 4852 : ABI_CHECK(rank /= self%nprocs, sjoin("Nobody has iq_ibz: ", itoa(iq_ibz)))
3492 : :
3493 : : ! Begin non-blocking communication for phfrq frequencies.
3494 [ + - + - : 36260 : if (self%my_rank == master) self%phfrq = self%phfreqs_qibz(:, iq_ibz)
- + - - +
+ ]
3495 : 4852 : call xmpi_ibcast(self%phfrq, master, self%comm, self%requests(1), ierr)
3496 : :
3497 : : ! Rotate eigvectors at q_ibz to get eigenvector at q_bz
3498 : : ! Don't test if umklapp == 0 because we use the periodic gauge:
3499 : : !
3500 : : ! phfreq(q+G) = phfreq(q) and eigvec(q) = eigvec(q+G)
3501 : : !
3502 [ + + + + ]: 4852 : isirr_q = isym_q == 1 .and. trev_q == 0
3503 : :
3504 [ + - ]: 4852 : if (self%my_rank == master) then
3505 : : ! I own the data --> operate on it
3506 [ + + ]: 4852 : if (isirr_q) then
3507 : : ! q in IBZ --> no rotation is needed.
3508 : : call phdispl_from_eigvec(cryst%natom, cryst%ntypat, cryst%typat, cryst%amu, &
3509 : 2014 : self%pheigvec_qibz(:,:,:,iq_ibz), self%displ_cart)
3510 : : else
3511 : : ! q in BZ --> rotate phonon eigenvectors.
3512 : : call pheigvec_rotate(cryst, self%qibz(:, iq_ibz), isym_q, trev_q, self%pheigvec_qibz(:,:,:,iq_ibz), &
3513 : 2838 : eigvec_qpt, self%displ_cart)
3514 : : end if
3515 : : end if
3516 : :
3517 : : ! Begin non-blocking bcast for displ_cart. Caller must wait (use phstore_wait)
3518 : 4852 : call xmpi_ibcast(self%displ_cart, master, self%comm, self%requests(2), ierr)
3519 : :
3520 : : end subroutine phstore_async_rotate
3521 : : !!***
3522 : :
3523 : : !----------------------------------------------------------------------
3524 : :
3525 : : !!****f* m_phonons/phstore_wait
3526 : : !! NAME
3527 : : !! phstore_wait
3528 : : !!
3529 : : !! FUNCTION
3530 : : !! Wait from non-blocking MPI BCAST started in phstore_async_rotate,
3531 : : !! returns phonon frequencies and displacements in Cartesian and reduced coordinates.
3532 : : !!
3533 : : !! INPUTS
3534 : :
3535 : 4852 : subroutine phstore_wait(self, cryst, phfrq, displ_cart, displ_red)
3536 : :
3537 : : !Arguments ------------------------------------
3538 : : class(phstore_t),intent(inout) :: self
3539 : : type(crystal_t),intent(in) :: cryst
3540 : : real(dp) ABI_ASYNC, intent(out) :: phfrq(self%natom3)
3541 : : real(dp) ABI_ASYNC, intent(out) :: displ_cart(2, 3, self%natom, self%natom3)
3542 : : real(dp),intent(out) :: displ_red(2, 3, self%natom, self%natom3)
3543 : :
3544 : : !Local variables ------------------------------
3545 : : integer :: ierr
3546 : : ! *************************************************************************
3547 : :
3548 [ + - ]: 4852 : if (.not. self%use_ifc_fourq) call xmpi_waitall(self%requests, ierr)
3549 [ + + ]: 31408 : phfrq = self%phfrq
3550 [ + + + + : 536968 : displ_cart = self%displ_cart
+ + + + ]
3551 : 4852 : call phdispl_cart2red(cryst%natom, cryst%gprimd, displ_cart, displ_red)
3552 : :
3553 : 4852 : end subroutine phstore_wait
3554 : : !!***
3555 : :
3556 : : !----------------------------------------------------------------------
3557 : :
3558 : : !!****f* m_phonons/test_phrotation
3559 : : !! NAME
3560 : : !! test_phrotation
3561 : : !!
3562 : : !! FUNCTION
3563 : : !! Test the symmetrization of the phonon eigenvalues and eigenvectors.
3564 : : !!
3565 : : !! INPUT
3566 : : !! cryst=Crystalline structure
3567 : : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
3568 : : !! ngqpt(3)=Divisions of the ab-initio q-mesh.
3569 : : !! qptopt=option for the generation of q points (defines whether spatial symmetries and/or time-reversal can be used)
3570 : : !! comm= MPI communicator
3571 : :
3572 : 0 : subroutine test_phrotation(ifc, cryst, qptopt, ngqpt, comm)
3573 : :
3574 : : !Arguments ------------------------------------
3575 : : type(ifc_type),intent(in) :: ifc
3576 : : type(crystal_t),intent(in) :: cryst
3577 : : integer,intent(in) :: qptopt, comm, ngqpt(3)
3578 : :
3579 : : !Local variables-------------------------------
3580 : : !scalars
3581 : : integer,parameter :: nqshft1 = 1, master = 0
3582 : : integer :: nqibz, iq_bz, iq_ibz, nqbz, ii, natom, natom3, ierr
3583 : : integer :: isym, itimrev, ierr_freq, ierr_eigvec, prtvol
3584 : : real(dp), parameter :: tol_phfreq_meV = tol3, tol_eigvec = tol6
3585 : : real(dp) :: maxerr_phfreq, err_phfreq, maxerr_eigvec ! err_eigvec
3586 : : logical :: isirr_q
3587 : : character(len=500) :: msg, fmt_freqs, fmt_eigvec
3588 : 0 : type(krank_t) :: qrank
3589 : : !arrays
3590 : : integer :: in_qptrlatt(3,3), new_qptrlatt(3,3), g0(3)
3591 : 0 : integer,allocatable :: bz2ibz(:,:), bz2ibz_listkk(:,:), toinv(:,:)
3592 : 0 : real(dp) :: qshift(3, nqshft1), phfrq(3*cryst%natom), work(3*cryst%natom)
3593 : 0 : real(dp) :: eigvec_out(2,3*cryst%natom,3*cryst%natom) !eigvec_ibz(2,3*cryst%natom,3*cryst%natom),
3594 : 0 : real(dp) :: eigvec_bz(2,3*cryst%natom,3*cryst%natom), displ_cart_qbz(2,3*cryst%natom,3*cryst%natom)
3595 : 0 : real(dp) :: d2cart(2,3*cryst%natom,3*cryst%natom), d2tmp(2,3*cryst%natom,3*cryst%natom)
3596 : 0 : real(dp),allocatable :: wtq_ibz(:), qbz(:,:), qibz(:,:), displ_cart(:,:,:,:),displ_red(:,:,:,:)
3597 : 0 : real(dp),allocatable :: phfreqs_qibz(:,:), displ_cart_ibz(:,:,:,:),eigvec_ibz(:,:,:,:)
3598 : : !************************************************************************
3599 : :
3600 [ # # ]: 0 : if (xmpi_comm_rank(comm) /= 0) return
3601 : :
3602 : 0 : call wrtout(std_out, sjoin(" Testing symmetrization of phonon frequencies and eigenvectors with ngqpt:", ltoa(ngqpt)), ch10)
3603 : 0 : prtvol = 1; natom = cryst%natom; natom3 = cryst%natom * 3
3604 : :
3605 : : ! Create a regular grid
3606 : 0 : in_qptrlatt = 0; in_qptrlatt(1, 1) = ngqpt(1); in_qptrlatt(2, 2) = ngqpt(2); in_qptrlatt(3, 3) = ngqpt(3)
3607 : 0 : qshift = zero
3608 : :
3609 : : call kpts_ibz_from_kptrlatt(cryst, in_qptrlatt, qptopt, nqshft1, qshift, &
3610 : 0 : nqibz, qibz, wtq_ibz, nqbz, qbz, new_kptrlatt=new_qptrlatt, bz2ibz=bz2ibz)
3611 [ # # ]: 0 : ABI_FREE(bz2ibz)
3612 : :
3613 : : !write(std_out, "(2(a, i0))")" nqibz: ", nqibz, ", nqbz:", nqbz
3614 : : !write(std_out, "(a)") " qibz_list:"
3615 : : !do iq_ibz=1,nqibz
3616 : : ! write(std_out, "(a)")trim(ltoa(qibz(:,iq_ibz)))
3617 : : !end do
3618 : : !call cryst%print(unit=std_out)
3619 : : !write(std_out, *)""
3620 : :
3621 : : ! Compute BZ --> IBZ mapping.
3622 [ # # # # : 0 : ABI_MALLOC(bz2ibz_listkk, (6, nqbz))
# # # # #
# # # ]
3623 : :
3624 : 0 : call qrank%from_kptrlatt(nqibz, qibz, in_qptrlatt, compute_invrank=.False.)
3625 : :
3626 [ # # ]: 0 : if (kpts_map("symrec", qptopt, cryst, qrank, nqbz, qbz, bz2ibz_listkk) /= 0) then
3627 : 0 : write(msg, '(3a)' ) "Error mapping BZ to IBZ",ch10,"The q-point could not be generated from a symmetrical one"
3628 : 0 : ABI_ERROR(msg)
3629 : : end if
3630 : :
3631 : 0 : call qrank%free()
3632 : :
3633 : : ! Compute ph freqs in the IBZ.
3634 [ # # # # : 0 : ABI_CALLOC(phfreqs_qibz, (natom3, nqibz))
# # # # #
# # # # #
# # # # #
# ]
3635 [ # # # # : 0 : ABI_CALLOC(displ_cart_ibz, (2, natom3, natom3, nqibz))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
3636 [ # # # # : 0 : ABI_CALLOC(eigvec_ibz, (2, natom3, natom3, nqibz))
# # # # #
# # # # #
# # # # #
# # # #
# ]
3637 : :
3638 [ # # ]: 0 : do iq_ibz=1,nqibz
3639 : : call ifc%fourq(cryst, qibz(:,iq_ibz), phfreqs_qibz(:,iq_ibz), displ_cart_ibz(:,:,:,iq_ibz), &
3640 : 0 : out_eigvec=eigvec_ibz(:,:,:,iq_ibz))
3641 : : end do
3642 : :
3643 [ # # # # : 0 : ABI_MALLOC(displ_cart, (2, 3, cryst%natom, natom3))
# # # # #
# # # # #
# # ]
3644 [ # # # # : 0 : ABI_MALLOC(displ_red, (2, 3, cryst%natom, natom3))
# # # # #
# # # ]
3645 : :
3646 [ # # # # : 0 : ABI_MALLOC(toinv, (4, cryst%nsym))
# # # # #
# # # ]
3647 : 0 : call sg_multable(cryst%nsym, cryst%symafm, cryst%symrel, ierr, toinv=toinv, tnons=cryst%tnons, tnons_tol=tol6)
3648 [ # # ]: 0 : ABI_CHECK(ierr == 0, "sg_multable returned ierr != 0")
3649 : :
3650 : : ! Precompute ph freqs in the BZ and compare with BZ
3651 : 0 : ierr_freq = 0; ierr_eigvec = 0
3652 [ # # ]: 0 : fmt_freqs = sjoin("(a, ", itoa(natom3), "(f7.3, 1x))")
3653 [ # # ]: 0 : fmt_eigvec = sjoin("(a, i0, 1x, a, ", itoa(natom3), "(f12.9, 1x))")
3654 : 0 : maxerr_phfreq = zero; maxerr_eigvec = zero
3655 : :
3656 [ # # ]: 0 : do iq_bz=1,nqbz
3657 : 0 : call ifc%fourq(cryst, qbz(:, iq_bz), phfrq, displ_cart, out_eigvec=eigvec_bz, out_d2cart=d2cart)
3658 : :
3659 : 0 : iq_ibz = bz2ibz_listkk(1, iq_bz); isym = bz2ibz_listkk(2, iq_bz)
3660 : 0 : itimrev = bz2ibz_listkk(6, iq_bz); g0 = bz2ibz_listkk(3:5, iq_bz)
3661 : : isirr_q = isym == 1 .and. itimrev == 0 .and. all(g0 == 0)
3662 : :
3663 : : ! Compare phfreqs within tol in meV.
3664 [ # # # # : 0 : err_phfreq = maxval(abs(phfrq - phfreqs_qibz(:, iq_ibz))) * Ha_meV
# # # # ]
3665 [ # # ]: 0 : if (err_phfreq > tol_phfreq_meV) then
3666 : 0 : maxerr_phfreq = max(maxerr_phfreq, err_phfreq)
3667 : 0 : write(std_out,*)" " // repeat("=", 92)
3668 [ # # # # ]: 0 : write(std_out, "(4a)")" qbz:", trim(ktoa(qbz(:, iq_bz))), " --> qibz:", trim(ktoa(qibz(:, iq_ibz)))
3669 [ # # ]: 0 : write(std_out, fmt_freqs)" w_bz :", phfrq * Ha_meV
3670 [ # # ]: 0 : write(std_out, fmt_freqs)" w_ibz:", phfreqs_qibz(:, iq_ibz) * Ha_meV
3671 : 0 : write(std_out,*)" err_phfreq (meV):", err_phfreq, " > tol: ", tol_phfreq_meV
3672 : 0 : write(std_out,*)" " // repeat("=", 92)
3673 : 0 : ierr_freq = ierr_freq + 1
3674 : : end if
3675 : :
3676 : : ! Rotate and compare eigenvectors
3677 : : call pheigvec_rotate(cryst, qibz(:, iq_ibz), isym, itimrev, eigvec_ibz(:,:,:,iq_ibz), &
3678 : 0 : eigvec_out, displ_cart_qbz)
3679 : :
3680 : : ! e^H D e = w**2 I
3681 : 0 : call massmult_and_breaksym(natom, cryst%ntypat, cryst%typat, cryst%amu, d2cart)
3682 : 0 : call cg_zgemm("N", "N", natom3, natom3, natom3, d2cart, eigvec_out, d2tmp)
3683 : 0 : call cg_zgemm("C", "N", natom3, natom3, natom3, eigvec_out, d2tmp, d2cart)
3684 [ # # ]: 0 : do ii=1,natom3
3685 : 0 : work(ii) = d2cart(1, ii, ii)
3686 : : end do
3687 [ # # ]: 0 : work = sqrt(abs(phfrq ** 2 - work)) * Ha_meV
3688 [ # # # # : 0 : if (maxval(work) > tol_phfreq_meV) ierr_eigvec = ierr_eigvec + 1
# # # # #
# ]
3689 [ # # # # : 0 : write(std_out, *) "max eig_diff [meV]: ", maxval(work)
# # # # ]
3690 : 0 : write(std_out, "(a)")" e^H D e (meV**2)"
3691 [ # # # # : 0 : d2cart = d2cart * Ha_meV ** 2
# # ]
3692 [ # # # # : 0 : call print_arr([std_out], reshape(cmplx(d2cart(1,:,:), d2cart(2,:,:), kind=dp), [natom3, natom3]))
# # # # #
# ]
3693 : :
3694 : : !err_eigvec = maxval(abs(eigvec_out - eigvec_bz))
3695 : : !if (err_eigvec > tol_eigvec) then
3696 : : ! maxerr_eigvec = max(maxerr_eigvec, err_eigvec)
3697 : : ! write(std_out, "(4a)")" qbz:", trim(ktoa(qbz(:, iq_bz))), " --> qibz: ", trim(ktoa(qibz(:, iq_ibz)))
3698 : : ! write(std_out, "(a,2(i0,1x),a)")" qbz image through isym, itimrev: ", isym, itimrev, trim(ltoa(cryst%tnons(:, isym)))
3699 : : ! write(std_out, "(a,l1)")" has_r0: ", any(cryst%indsym(1:3, isym, :) /= 0)
3700 : : ! write(std_out, *) "err_eigvec ", err_eigvec, " > tol:", tol_eigvec
3701 : : ! do ii=1,natom3
3702 : : ! if (all(abs(eigvec_out(:,:,ii) - eigvec_bz(:,:,ii)) < tol_eigvec)) cycle
3703 : : ! if (prtvol > 0) then
3704 : : ! write(std_out, "(a, 2(f12.9,1x))") " diff (re/im): ", &
3705 : : ! maxval(abs(eigvec_out(1,:,ii) - eigvec_bz(1,:,ii))), maxval(abs(eigvec_out(2,:,ii) - eigvec_bz(2,:,ii)))
3706 : : ! write(std_out, fmt_eigvec) " mode: ", ii, "re_sym: ", eigvec_out(1,:,ii)
3707 : : ! write(std_out, fmt_eigvec) " mode: ", ii, "re_bz : ", eigvec_bz(1,:,ii)
3708 : : ! write(std_out, fmt_eigvec) " mode: ", ii, "im_sym: ", eigvec_out(2,:,ii)
3709 : : ! write(std_out, fmt_eigvec) " mode: ", ii, "im_bz : ", eigvec_bz(2,:,ii)
3710 : : ! end if
3711 : : ! end do
3712 : : ! write(std_out,*)" " // repeat("=", 92)
3713 : : ! ierr_eigvec = ierr_eigvec + 1
3714 : : !end if
3715 : :
3716 : : end do
3717 : :
3718 : 0 : write(std_out,*)" === Final results ==="
3719 : 0 : write(std_out,*)" maxerr_phfreq (meV): ", maxerr_phfreq
3720 : 0 : write(std_out,*)" percentage of erroneous q-points for phfreq: ", ierr_freq / (one * nqbz) * 100, "%"
3721 : 0 : write(std_out,*)""
3722 : : !write(std_out,*)" maxerr_eigvec: ", maxerr_eigvec
3723 : 0 : write(std_out,*)" percentage of q-points for eigvec: ", ierr_eigvec / (one * nqbz) * 100, "%"
3724 : 0 : write(std_out,*)""
3725 : :
3726 [ # # ]: 0 : if (ierr_freq /= 0) then
3727 : 0 : ABI_ERROR("Wrong symmetrization in phonon eigenvalues.")
3728 [ # # ]: 0 : else if (ierr_eigvec /= 0) then
3729 : 0 : ABI_ERROR("Wrong symmetrization in phonon eigenvectors.")
3730 : : else
3731 : 0 : write(std_out, "(a)")" ALL OK: No error detected!"
3732 : : end if
3733 : :
3734 [ # # ]: 0 : ABI_FREE(displ_cart)
3735 [ # # ]: 0 : ABI_FREE(displ_red)
3736 [ # # ]: 0 : ABI_SFREE(bz2ibz_listkk)
3737 [ # # ]: 0 : ABI_SFREE(qbz)
3738 [ # # ]: 0 : ABI_SFREE(qibz)
3739 [ # # ]: 0 : ABI_SFREE(wtq_ibz)
3740 [ # # ]: 0 : ABI_FREE(phfreqs_qibz)
3741 [ # # ]: 0 : ABI_FREE(displ_cart_ibz)
3742 [ # # ]: 0 : ABI_FREE(eigvec_ibz)
3743 [ # # ]: 0 : ABI_FREE(toinv)
3744 : :
3745 [ # # # # : 0 : end subroutine test_phrotation
# # # # #
# # # #
# ]
3746 : : !!***
3747 : :
3748 [ - + - - : 1253 : end module m_phonons
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ - + - +
- + ]
3749 : : !!***
|
