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