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 77 : 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 77 : unitname='Ha'
297 :
298 129 : select case (PHdos%prtdos)
299 : case (1)
300 : write(msg_method,'(a,es16.8,2a,i0)')&
301 52 : '# 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 77 : ABI_ERROR(sjoin(" Wrong prtdos: ",itoa(PHdos%prtdos)))
306 : end select
307 :
308 : ! Open external file and write results
309 77 : if (open_file(fname,msg,newunit=unt,form="formatted",action="write") /= 0) then
310 0 : ABI_ERROR(msg)
311 : end if
312 77 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom type projected DOS'
313 77 : call wrtout(unt,msg)
314 77 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
315 77 : call wrtout(unt,msg)
316 77 : call wrtout(unt,msg_method)
317 77 : write(msg,'(5a)')'# ',ch10,'# omega PHDOS INT_PHDOS PJDOS[atom_type=1] INT_PJDOS[atom_type=1] ... ',ch10,'# '
318 77 : call wrtout(unt,msg)
319 144574 : do io=1,PHdos%nomega
320 144497 : write(unt,'(3es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io),PHdos%phdos_int(io)
321 345113 : do itype=1,PHdos%ntypat
322 345113 : write(unt,'(2es17.8,2x)',advance='NO')PHdos%pjdos_type(io,itype),PHdos%pjdos_type_int(io,itype)
323 : end do
324 144574 : write(unt,*)
325 : end do
326 77 : close(unt)
327 :
328 77 : fname_by_atom = trim(fname) // "_by_atom"
329 77 : 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 77 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states and atom projected DOS'
333 77 : call wrtout(unt_by_atom,msg)
334 77 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in states/',unitname
335 77 : call wrtout(unt_by_atom,msg)
336 77 : call wrtout(unt_by_atom,msg_method)
337 77 : write(msg,'(5a)')'# ',ch10,'# omega PHDOS PJDOS[atom=1] PJDOS[atom=2] ... ',ch10,'# '
338 77 : call wrtout(unt_by_atom,msg)
339 144574 : do io=1,PHdos%nomega
340 144497 : write(unt_by_atom,'(2es17.8)',advance='NO')PHdos%omega(io),PHdos%phdos(io)
341 498707 : do iatom=1,PHdos%natom
342 1561337 : write(unt_by_atom,'(1es17.8,2x)',advance='NO') sum(PHdos%pjdos(io,1:3,iatom))
343 : end do
344 144574 : write(unt_by_atom,*)
345 : end do
346 77 : close(unt_by_atom)
347 :
348 77 : fname_msqd = trim(fname) // "_msqd"
349 77 : if (open_file(fname_msqd,msg,newunit=unt_msqd,form="formatted",action="write") /= 0) then
350 0 : ABI_ERROR(msg)
351 : end if
352 77 : write(msg,'(3a)')'# ',ch10,'# Phonon density of states weighted msq displacement matrix (set to zero below 1e-12)'
353 77 : call wrtout(unt_msqd,msg)
354 77 : write(msg,'(6a)')'# ',ch10,'# Energy in ',unitname,', DOS in bohr^2 states/',unitname
355 77 : call wrtout(unt_msqd,msg)
356 77 : call wrtout(unt_msqd,msg_method)
357 77 : write(msg,'(5a)')'# ',ch10,'# omega MSQDisp[atom=1, xx, yy, zz, yz, xz, xy] MSQDisp[atom=2, xx, yy,...] ... ',ch10,'# '
358 77 : call wrtout(unt_msqd,msg)
359 144574 : do io=1,PHdos%nomega
360 144497 : write(unt_msqd,'(2es17.8)',advance='NO')PHdos%omega(io)
361 498707 : do iatom=1,PHdos%natom
362 4604730 : tens = PHdos%msqd_dos_atom(io,:,:,iatom)
363 4604730 : where (abs(tens) < tol12)
364 : tens = zero
365 : end where
366 : write(unt_msqd,'(6es17.8,2x)',advance='NO') &
367 354210 : tens(1,1), tens(2,2), tens(3,3), &
368 852917 : tens(2,3), tens(1,3), tens(1,2)
369 : end do
370 144574 : write(unt_msqd,*)
371 : end do
372 77 : close(unt_msqd)
373 :
374 77 : 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 173 : subroutine phdos_free(PHdos)
600 :
601 : !Arguments -------------------------------
602 : class(phdos_t),intent(inout) ::PHdos
603 : ! *************************************************************************
604 :
605 173 : ABI_SFREE(PHdos%atom_mass)
606 173 : ABI_SFREE(PHdos%normal_vec_dmm)
607 173 : ABI_SFREE(PHdos%omega)
608 173 : ABI_SFREE(PHdos%phdos)
609 173 : ABI_SFREE(PHdos%phdos_int)
610 173 : ABI_SFREE(PHdos%phdos_dmm)
611 173 : ABI_SFREE(PHdos%pjdos)
612 173 : ABI_SFREE(PHdos%pjdos_int)
613 173 : ABI_SFREE(PHdos%pjdos_type)
614 173 : ABI_SFREE(PHdos%pjdos_type_int)
615 173 : ABI_SFREE(PHdos%pjdos_rc_type)
616 173 : ABI_SFREE(PHdos%msqd_dos_atom)
617 :
618 173 : 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 173 : 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 173 : phdos%ntypat = crystal%ntypat
652 173 : phdos%natom = crystal%natom
653 173 : phdos%prtdos = prtdos
654 173 : phdos%dossmear = dossmear
655 173 : phdos%omega_step = dosdeltae
656 173 : phdos%n_normal_vec_dmm = 7
657 173 : 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 173 : if (wminmax(2) > wminmax(1)) then
662 66 : phdos%omega_min = wminmax(1)
663 66 : phdos%omega_max = wminmax(2)
664 : else
665 107 : phdos%omega_min = ifc%omega_minmax(1)
666 107 : phdos%omega_max = ifc%omega_minmax(2)
667 : end if
668 : ! Must be consistent with mesh computed in tetra routines!
669 173 : phdos%nomega = nint((phdos%omega_max - phdos%omega_min) / phdos%omega_step) + 1
670 : ! Ensure Simpson integration will be ok
671 173 : phdos%nomega = max(6, phdos%nomega)
672 :
673 : ! Build frequency mesh.
674 519 : ABI_MALLOC(phdos%omega, (phdos%nomega))
675 194872 : do io=1,phdos%nomega
676 194872 : phdos%omega(io) = phdos%omega_min + phdos%omega_step * (io - 1)
677 : end do
678 173 : phdos%omega_min = phdos%omega(1)
679 173 : phdos%omega_max = phdos%omega(phdos%nomega)
680 :
681 : ! Allocate arrays that depend on nomega and set them to zero.
682 5363 : ABI_CALLOC(phdos%normal_vec_dmm, (3, phdos%n_normal_vec_dmm))
683 195045 : ABI_CALLOC(phdos%phdos, (phdos%nomega))
684 195045 : ABI_CALLOC(phdos%phdos_int, (phdos%nomega))
685 1364796 : ABI_CALLOC(phdos%phdos_dmm, (phdos%nomega,phdos%n_normal_vec_dmm))
686 1369680 : ABI_CALLOC(phdos%pjdos, (phdos%nomega, 3, crystal%natom))
687 1369507 : ABI_CALLOC(phdos%pjdos_int, (phdos%nomega, 3, crystal%natom))
688 4108199 : ABI_CALLOC(phdos%msqd_dos_atom, (phdos%nomega, 3, 3, crystal%natom))
689 519 : ABI_MALLOC(phdos%atom_mass, (crystal%natom))
690 1086 : phdos%atom_mass = crystal%amu(crystal%typat(:)) * amu_emass
691 :
692 173 : 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 173 : 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 173 : 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 173 : integer,allocatable :: bz2ibz_smap(:,:), bz2ibz(:)
771 : real(dp) :: veloc(3), speedofsound(3),speedofsound_(3)
772 346 : real(dp) :: displ(2*3*Crystal%natom*3*Crystal%natom)
773 519 : real(dp) :: eigvec(2,3,Crystal%natom,3*Crystal%natom),phfrq(3*Crystal%natom),phangmom(3,3*Crystal%natom)
774 346 : 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 346 : real(dp) :: symcart(3,3,crystal%nsym), syme2_xyza(3, crystal%natom)
777 173 : real(dp),allocatable :: full_eigvec(:,:,:,:,:),full_phfrq(:,:),full_phangmom(:,:,:),new_shiftq(:,:), full_veloc(:,:,:)
778 173 : real(dp),allocatable :: qbz(:,:),qibz(:,:),tmp_phfrq(:) !, work_msqd(:,:,:,:)
779 173 : real(dp),allocatable :: wtq_ibz(:),xvals(:), gvals_wtq(:), wdt(:,:), energies(:)
780 : ! *********************************************************************
781 :
782 : DBG_ENTER("COLL")
783 :
784 173 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
785 519 : units = [std_out, ab_out]
786 173 : prtout_ = .true.
787 173 : if (present(prtout)) then
788 88 : prtout_ = prtout
789 : end if
790 :
791 : ! Consistency check.
792 173 : 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 173 : dosdeltae = dosdeltae_in; refine_dosdeltae = .false.
796 173 : 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 173 : 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 173 : call cwtime(cpu_all, wall_all, gflops_all, "start")
805 :
806 : ! Get symmetries in cartesian coordinates
807 7167 : do isym=1,crystal%nsym
808 7167 : call symredcart(crystal%rprimd,crystal%gprimd,symcart(:,:,isym),crystal%symrel(:,:,isym))
809 : end do
810 :
811 173 : natom = crystal%natom
812 173 : call phdos_malloc(phdos, crystal, ifc, dosdeltae, dossmear, wminmax, prtdos)
813 173 : nomega = phdos%nomega
814 :
815 4613 : veloc_1q = zero
816 : ! TODO: add input variable and adapt these directions to be read in from input file
817 173 : if (phdos%n_normal_vec_dmm >= 7) then
818 692 : phdos%normal_vec_dmm(:,1) = [one, zero, zero]
819 692 : phdos%normal_vec_dmm(:,2) = [zero, one, zero]
820 692 : phdos%normal_vec_dmm(:,3) = [zero, zero, one]
821 692 : phdos%normal_vec_dmm(:,4) = [one, one, zero]
822 692 : phdos%normal_vec_dmm(:,5) = [one, zero, one]
823 692 : phdos%normal_vec_dmm(:,6) = [zero, one, one]
824 692 : phdos%normal_vec_dmm(:,7) = [one, one, one]
825 : end if
826 519 : ABI_MALLOC(gvals_wtq, (nomega))
827 346 : ABI_MALLOC(xvals, (nomega))
828 :
829 : ! Parameters defining the gaussian approximant.
830 173 : if (prtdos == 1) then
831 : ! TODO: use gaussian and update reference files.
832 97 : gaussprefactor = one / (dossmear * sqrt(two_pi))
833 97 : gaussfactor = one / (sqrt2 * dossmear)
834 97 : write(msg, '(4a,f8.5,2a,f8.5,a,i0)') ch10, &
835 97 : ' phdos_init: calculating phonon DOS using gaussian method:', ch10, &
836 97 : ' gaussian smearing [meV] = ', dossmear * Ha_meV, ch10, &
837 194 : ' 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 173 : if (prtout_) call wrtout(std_out, msg)
844 :
845 : ! This call will set %nqibz and IBZ and BZ arrays
846 173 : 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 173 : 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 173 : 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 2249 : phdos%qptrlatt = new_qptrlatt
856 692 : phdos%shiftq(:) = new_shiftq(:, 1) ! only one shift in output
857 :
858 173 : if (my_rank == master) then
859 173 : write(msg, "(3a, i0)")" DOS ngqpt: ", trim(ltoa(dos_ngqpt)), ", qptopt: ", my_qptopt
860 173 : if (prtout_) call wrtout(std_out, msg)
861 173 : write(msg, "(2(a, i0))")" Number of q-points in the IBZ: ", phdos%nqibz, ", number of MPI processes: ", nprocs
862 173 : if (prtout_) call wrtout(std_out, msg)
863 : end if
864 : !call cwtime_report(" kpts_ibz_from_kptrlatt", cpu, wall, gflops)
865 :
866 173 : 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 173 : ABI_FREE(bz2ibz_smap)
889 173 : 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 173 : nsmallq = zero; speedofsound = zero
898 173 : wminmax = [huge(one), -huge(one)]; count_wminmax = 0
899 :
900 227454 : do iq_ibz=1,phdos%nqibz
901 227281 : 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 227281 : if (phdos%n_normal_vec_dmm > 0) then
907 227281 : 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 1862849 : wminmax(1) = min(wminmax(1), minval(phfrq))
913 227281 : if (wminmax(1) < phdos%omega(1)) count_wminmax(1) = count_wminmax(1) + 1
914 1862849 : wminmax(2) = max(wminmax(2), maxval(phfrq))
915 227281 : if (wminmax(2) > phdos%omega(nomega)) count_wminmax(2) = count_wminmax(2) + 1
916 :
917 909124 : normq = sum(qibz(:,iq_ibz) ** 2)
918 227281 : if (normq < max_smallq .and. normq > tol6) then
919 17137 : call phdos_calc_vsound(eigvec, crystal%gmet, natom, phfrq, qibz(:,iq_ibz), speedofsound_)
920 68548 : speedofsound = speedofsound + speedofsound_ * wtq_ibz(iq_ibz)
921 17137 : nsmallq = nsmallq + wtq_ibz(iq_ibz)
922 : end if
923 :
924 227281 : dos_maxmode_ = 3*natom
925 227281 : if (present(dos_maxmode)) then
926 1982 : if (dos_maxmode > 0 .and. dos_maxmode < 3*natom) then
927 227281 : dos_maxmode_ = dos_maxmode
928 : end if
929 : end if
930 :
931 173 : select case (prtdos)
932 : case (1)
933 1531548 : do imode=1,dos_maxmode_
934 : ! Precompute \delta(w - w_{qnu}) * weight(q)
935 1085258169 : xvals = (phdos%omega(:) - phfrq(imode)) * gaussfactor
936 1083938277 : 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 1083938277 : 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 18210180 : syme2_xyza = zero
949 5542464 : do iat=1,natom
950 72584064 : do isym=1,my_nsym
951 67041600 : jat = crystal%indsym(4,isym,iat)
952 : syme2_xyza(:,jat) = syme2_xyza(:,jat) + &
953 : matmul(symcart(:,:,isym), eigvec(1,:,iat,imode)) ** 2 + &
954 1881387372 : matmul(symcart(:,:,isym), eigvec(2,:,iat,imode)) ** 2
955 : end do
956 : end do
957 18210180 : syme2_xyza = syme2_xyza / my_nsym
958 :
959 : ! Accumulate PJDOS
960 5542464 : do iat=1,natom
961 18210180 : do idir=1,3
962 10147882227 : 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 5542464 : do iat=1,natom
971 4222572 : msqd_atom_tmp = zero
972 16890288 : do idir=1,3
973 54893436 : 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 50670864 : + 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 72584064 : do isym=1, my_nsym
986 : !temp_33 = matmul( (symcart(:,:,isym)), matmul(msqd_atom_tmp, transpose(symcart(:,:,isym))) )
987 : ! MG Version
988 6100785600 : temp_33 = matmul( (transpose(symcart(:,:,isym))), matmul(msqd_atom_tmp, symcart(:,:,isym)) )
989 871540800 : temp_33 = temp_33 / my_nsym
990 67041600 : jat = crystal%indsym(4,isym,iat)
991 272388972 : do idir=1,3
992 871540800 : do jdir=1,3
993 : phdos%msqd_dos_atom(:,idir,jdir,jat) = phdos%msqd_dos_atom(:,idir,jdir,jat) + &
994 >57027*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 5279568 : veloc = veloc_1q(:,imode)
1003 10770792 : 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 199451532 : do isym=1,my_nsym
1010 3052635852 : projfact = projfact + abs(dot_product(matmul(symcart(:,:,isym),veloc),phdos%normal_vec_dmm(:,icart)))
1011 : end do
1012 7588887831 : 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 227281 : ABI_ERROR(sjoin("Wrong value for prtdos:", itoa(prtdos)))
1027 : end select
1028 : end do ! iq_ibz
1029 :
1030 173 : ABI_FREE(qbz)
1031 173 : ABI_FREE(gvals_wtq)
1032 173 : ABI_FREE(xvals)
1033 :
1034 173 : call xmpi_sum_master(nsmallq, master, comm, ierr)
1035 173 : call xmpi_sum_master(speedofsound, master, comm, ierr)
1036 :
1037 : !call cwtime_report(" phdos", cpu, wall, gflops)
1038 :
1039 173 : if (my_rank == master) then
1040 173 : if (nsmallq > tol10) then
1041 : ! Write info about speed of sound
1042 684 : speedofsound = speedofsound / nsmallq
1043 : write (msg,'(a,E20.10,3a,F16.4,2a)') &
1044 684 : ' Average speed of sound partial sums: ', third*sum(speedofsound), ' (at units)',ch10, &
1045 855 : '- = ', third*sum(speedofsound) * Bohr_Ang * 1.d-13 / Time_Sec, ' [km/s]',ch10
1046 171 : if (prtout_) call wrtout(units, msg)
1047 :
1048 : ! Debye frequency = vs * (6 pi^2 natom / ucvol)**1/3
1049 684 : debyefreq = third*sum(speedofsound) * (six*pi**2/crystal%ucvol)**(1./3.)
1050 : write (msg,'(a,E20.10,3a,E20.10,a)') &
1051 171 : ' Debye frequency from partial sums: ', debyefreq, ' (Ha)',ch10, &
1052 342 : '- = ', debyefreq*Ha_THz, ' (THz)'
1053 171 : if (prtout_) call wrtout(units, msg)
1054 :
1055 : ! Debye temperature = hbar * Debye frequency / kb
1056 171 : write (msg,'(a,E20.10,2a)') '-Debye temperature from partial sums: ', debyefreq*Ha_K, ' (K)', ch10
1057 171 : if (prtout_) call wrtout(units, msg)
1058 : end if
1059 : end if
1060 :
1061 173 : 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 97 : 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 173 : call xmpi_sum(count_wminmax, comm, ierr)
1210 346 : call xmpi_min(wminmax(1), rtmp, comm, ierr); wminmax(1) = rtmp
1211 346 : call xmpi_max(wminmax(2), rtmp, comm, ierr); wminmax(2) = rtmp
1212 :
1213 173 : call xmpi_sum(phdos%phdos, comm, ierr)
1214 173 : call xmpi_sum(phdos%phdos_dmm, comm, ierr)
1215 173 : call xmpi_sum(phdos%msqd_dos_atom, comm, ierr)
1216 173 : call xmpi_sum(phdos%pjdos, comm, ierr)
1217 173 : 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 1364277 : phdos%phdos_dmm = phdos%phdos_dmm * two_pi**3 / crystal%ucvol
1225 :
1226 : ! ===============================
1227 : ! === Compute Integrated PDOS ===
1228 : ! ===============================
1229 862333 : ABI_CALLOC(phdos%pjdos_rc_type, (nomega, 3, crystal%ntypat))
1230 287825 : ABI_CALLOC(phdos%pjdos_type, (nomega, crystal%ntypat))
1231 287652 : ABI_CALLOC(phdos%pjdos_type_int, (nomega, crystal%ntypat))
1232 :
1233 543 : do iat=1,natom
1234 370 : itype = crystal%typat(iat)
1235 456206 : do io=1,phdos%nomega
1236 1823344 : phdos%pjdos_rc_type(io,:,itype) = phdos%pjdos_rc_type(io,:,itype) + phdos%pjdos(io,:,iat)
1237 1823714 : phdos%pjdos_type(io,itype) = phdos%pjdos_type(io,itype) + sum(phdos%pjdos(io,:,iat))
1238 : end do
1239 543 : 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 173 : if (prtdos == 1) then
1249 97 : call simpson_int(phdos%nomega,phdos%omega_step, phdos%phdos, phdos%phdos_int)
1250 318 : do iat=1,natom
1251 981 : do idir=1,3
1252 884 : call simpson_int(phdos%nomega, phdos%omega_step, phdos%pjdos(:,idir,iat), phdos%pjdos_int(:,idir,iat))
1253 : end do
1254 : end do
1255 216 : do itype=1,crystal%ntypat
1256 216 : 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 173 : ABI_FREE(qibz)
1261 173 : ABI_FREE(wtq_ibz)
1262 :
1263 173 : call cwtime_report(" phdos_init", cpu_all, wall_all, gflops_all)
1264 : DBG_EXIT("COLL")
1265 :
1266 519 : 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 17374 : 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 17374 : imode_acoustic = 0
2140 :
2141 123013 : do imode = 1, 3*natom
2142 : ! Check if this mode is acoustic like: scalar product of all displacement vectors are collinear
2143 423768 : isacoustic = 1
2144 : ! Find reference atom with non-zero displacement
2145 423768 : do iatom=1,natom
2146 3286929 : 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 423768 : do iatom = 1, natom
2150 3181290 : if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, imode)&
2151 211289 : *eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, imode)) < tol16 ) isacoustic = 0
2152 : end do
2153 105639 : if (isacoustic == 0) cycle
2154 52203 : 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 678639 : qtmp = matmul(gmet, qphon)
2159 208812 : qnormcart = two * pi * sqrt(sum(qphon*qtmp))
2160 123013 : speedofsound(imode_acoustic) = phfrq(imode) / qnormcart
2161 : end do
2162 :
2163 17374 : 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 1229 : end module m_phonons
3727 : !!***
|
