Line data Source code
1 : !!****m* ABINIT/m_gruneisen
2 : !! NAME
3 : !! m_gruneisen
4 : !!
5 : !! FUNCTION
6 : !! Objects and methods to compute Gruneisen parameters with central finite differences
7 : !! of dynamical matrices obtained with different unit cell volumes.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2011-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_gruneisen
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_crystal
30 : use m_htetra
31 : use m_ddb
32 : use m_ddb_hdr
33 : use m_ifc
34 : use m_cgtools
35 : use m_nctk
36 : use netcdf
37 :
38 : use m_io_tools, only : get_unit, open_file
39 : use m_time, only : cwtime, cwtime_report
40 : use m_fstrings, only : sjoin, itoa, ltoa, ftoa, strcat
41 : use m_numeric_tools, only : central_finite_diff, arth
42 : use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt
43 : use m_bz_mesh, only : kpath_t
44 : use m_anaddb_dataset, only : anaddb_dataset_type
45 : use m_dynmat, only : massmult_and_breaksym, dfpt_phfrq, gtdyn9
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : !!****t* m_gruneisen/gruns_t
53 : !! NAME
54 : !! gruns_t
55 : !!
56 : !! FUNCTION
57 : !! Contains the interatomic force constants for different volumes.
58 : !! Provides methods to compute phonon bandstructures and gruneisen parameters.
59 : !!
60 : !! SOURCE
61 :
62 : type,public :: gruns_t
63 :
64 : integer :: natom3
65 : ! 3 * natom
66 :
67 : integer :: nvols
68 : ! Number of volumes.
69 :
70 : integer :: iv0
71 : ! Index of the DDB file corresponding to the equilibrium volume V0.
72 :
73 : real(dp) :: v0
74 : ! Equilibrium volume.
75 :
76 : real(dp) :: delta_vol
77 : ! Uniform grid spacing for finite difference.
78 :
79 : type(crystal_t),allocatable :: cryst_vol(:)
80 : ! cryst_vol(nvols)
81 : ! crystalline structure for the different volumes.
82 :
83 : type(ddb_type),allocatable :: ddb_vol(:)
84 : ! dbb_vol(nvols)
85 : ! DDB objects for the different volumes.
86 :
87 : type(ifc_type),allocatable :: ifc_vol(:)
88 : ! ifc_vol(nvols)
89 : ! interatomic force constants for the different volumes.
90 :
91 : end type gruns_t
92 :
93 : public :: gruns_new ! Constructor.
94 : public :: gruns_qpath ! Compute Gruneisen parameters on a q-path.
95 : public :: gruns_qmesh ! Compute Gruneisen parameters on a q-mesh.
96 : public :: gruns_free ! Release memory.
97 : public :: gruns_anaddb ! Driver routine called in anaddb.
98 : !!***
99 :
100 : contains !===========================================================
101 : !!***
102 :
103 : !----------------------------------------------------------------------
104 :
105 : !!****f* m_gruneisen/gruns_new
106 : !! NAME
107 : !! gruns_new
108 : !!
109 : !! FUNCTION
110 : !! Construct new object from a list of DDB files.
111 : !!
112 : !! INPUTS
113 : !! ddb_filepaths(:)=Paths of the DDB files (must be ordered by volume)
114 : !! inp<anaddb_dataset_type>=Anaddb dataset with input variables
115 : !! comm=MPI communicator
116 : !!
117 : !! SOURCE
118 :
119 1 : type(gruns_t) function gruns_new(ddb_filepaths, inp, comm) result(new)
120 :
121 : !Arguments ------------------------------------
122 : integer,intent(in) :: comm
123 : type(anaddb_dataset_type),intent(in) :: inp
124 : !arrays
125 : character(len=*),intent(in) :: ddb_filepaths(:)
126 :
127 : !Local variables-------------------------------
128 : integer,parameter :: master=0
129 : integer :: ivol,iblock,natom,ddbun, nprocs,my_rank,ierr
130 : character(len=500) :: msg
131 52 : type(ddb_hdr_type) :: ddb_hdr
132 : !arrays
133 : real(dp) :: dielt(3,3)
134 1 : real(dp),allocatable :: zeff(:,:,:), qdrp_cart(:,:,:,:)
135 :
136 : ! ************************************************************************
137 :
138 1 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
139 :
140 1 : new%nvols = size(ddb_filepaths)
141 57 : ABI_MALLOC(new%cryst_vol, (new%nvols))
142 6 : ABI_MALLOC(new%ddb_vol, (new%nvols))
143 11 : ABI_MALLOC(new%ifc_vol, (new%nvols))
144 :
145 1 : call wrtout(ab_out, "Computation of Gruneisen parameter with central finite difference:")
146 :
147 1 : ddbun = get_unit()
148 4 : do ivol=1,new%nvols
149 3 : call wrtout(ab_out, sjoin(" Reading DDB file:", ddb_filepaths(ivol)))
150 :
151 3 : call new%ddb_vol(ivol)%from_file(ddb_filepaths(ivol), ddb_hdr, new%cryst_vol(ivol), comm)
152 3 : call new%ddb_vol(ivol)%set_brav(inp%brav)
153 3 : natom = ddb_hdr%natom
154 3 : call ddb_hdr%free()
155 :
156 3 : if (my_rank == master) then
157 3 : call new%cryst_vol(ivol)%print(header=sjoin("Structure for ivol:", itoa(ivol)), unit=ab_out, prtvol=-1)
158 : end if
159 :
160 : ! Get Dielectric Tensor and Effective Charges
161 : ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
162 9 : ABI_MALLOC(zeff, (3,3,natom))
163 369 : ABI_CALLOC(qdrp_cart, (3,3,3,natom))
164 3 : iblock = new%ddb_vol(ivol)%get_dielt_zeff(new%cryst_vol(ivol), inp%rfmeth, inp%chneut, inp%selectz, dielt, zeff)
165 3 : if (iblock == 0) then
166 0 : call wrtout(ab_out, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_filepaths(ivol)))
167 0 : call wrtout(ab_out, "Values initialized with zeros")
168 : else
169 3 : call wrtout(ab_out, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_filepaths(ivol)))
170 : end if
171 :
172 : call new%ifc_vol(ivol)%init(new%cryst_vol(ivol), new%ddb_vol(ivol),&
173 : inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,zeff,&
174 3 : qdrp_cart,inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,inp%sys_dim,comm)
175 3 : ABI_FREE(zeff)
176 4 : ABI_FREE(qdrp_cart)
177 : end do
178 :
179 : ! Consistency check
180 : ! TODO: Add more tests.
181 1 : ABI_CHECK(any(new%nvols == [3, 5, 7, 9]), "Central finite difference requires [3,5,7,9] DDB files")
182 :
183 1 : new%natom3 = 3 * new%cryst_vol(1)%natom
184 1 : new%iv0 = 1 + new%nvols / 2
185 1 : new%v0 = new%cryst_vol(new%iv0)%ucvol
186 1 : new%delta_vol = new%cryst_vol(new%iv0+1)%ucvol - new%cryst_vol(new%iv0)%ucvol
187 :
188 1 : ierr = 0
189 4 : do ivol=1,new%nvols
190 4 : if (abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol) > tol4) then
191 0 : write(std_out,*)"ucvol, delta_vol, diff_vol", new%cryst_vol(ivol)%ucvol, new%delta_vol, &
192 0 : abs(new%cryst_vol(1)%ucvol + new%delta_vol * (ivol-1) - new%cryst_vol(ivol)%ucvol)
193 0 : ierr = ierr + 1
194 : end if
195 : end do
196 1 : if (ierr /= 0) then
197 0 : msg = ltoa([(new%cryst_vol(ivol)%ucvol, ivol=1,new%nvols)])
198 0 : ABI_ERROR(sjoin("Gruneisen calculations requires linear mesh of volumes but received:", msg))
199 : end if
200 :
201 2 : end function gruns_new
202 : !!***
203 :
204 : !----------------------------------------------------------------------
205 :
206 : !!****f* m_gruneisen/gruns_fourq
207 : !! NAME
208 : !! gruns_fourq
209 : !!
210 : !! FUNCTION
211 : !! Compute gruneisen parameters at an arbitrary q-point.
212 : !!
213 : !! INPUTS
214 : !! qpt(3)=q-point in reduced coordinates.
215 : !!
216 : !! OUTPUT
217 : !! wvols(3*natom, nvols) = Phonon frequencies for the different volumen.
218 : !! gvals(3*natom)=Gruneisen parameters evaluated at V0.
219 : !! dwdq(3,3*natom)=Group velocities at V0 in Cartesian coordinates.
220 : !! phdispl_cart(2, natom3, natom3, nvols)=Phonon displacement in Cartesian coordinates for the different volumes
221 : !!
222 : !! NOTES
223 : !!
224 : !! The Gruneisen parameter is given by:
225 : !!
226 : !! gamma(q,nu) = - (V / w(q,nu)) dw(q,nu)/dV
227 : !!
228 : !! Using w*2 = <u|D|u> and the Hellmann-Feynmann theorem, one obtains:
229 : !!
230 : !! gamma(q,nu) = - (V / 2 w(q,nu)**2) <u(q,nu)|dD(q)/dV|u(q,nu)>
231 : !!
232 : !! The derivative dD/dV is computed via central finite difference.
233 : !!
234 : !! SOURCE
235 :
236 228 : subroutine gruns_fourq(gruns, qpt, wvols, gvals, dwdq, phdispl_cart)
237 :
238 : !Arguments ------------------------------------
239 : !scalars
240 : class(gruns_t),intent(in) :: gruns
241 : !arrays
242 : real(dp),intent(in) :: qpt(3)
243 : real(dp),intent(out) :: wvols(gruns%natom3,gruns%nvols),gvals(gruns%natom3),dwdq(3,gruns%natom3)
244 : real(dp),intent(out) :: phdispl_cart(2, gruns%natom3, gruns%natom3, gruns%nvols)
245 :
246 : !Local variables-------------------------------
247 : !scalars
248 : integer :: ivol,natom3,nu
249 : real(dp) :: fact
250 : !arrays
251 : real(dp) :: dot(2)
252 456 : real(dp) :: eigvec(2,gruns%natom3,gruns%natom3,gruns%nvols),d2cart(2,gruns%natom3,gruns%natom3,gruns%nvols)
253 456 : real(dp) :: dddv(2,gruns%natom3,gruns%natom3)
254 456 : real(dp) :: omat(2,gruns%natom3, gruns%natom3)
255 :
256 : ! ************************************************************************
257 :
258 228 : natom3 = gruns%natom3
259 912 : do ivol=1,gruns%nvols
260 684 : if (ivol == gruns%iv0) then
261 : ! Compute group velocities for V=V0
262 : call gruns%ifc_vol(ivol)%fourq(gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), &
263 228 : out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol), dwdq=dwdq)
264 : else
265 : call gruns%ifc_vol(ivol)%fourq(gruns%cryst_vol(ivol), qpt, wvols(:,ivol), phdispl_cart(:,:,:,ivol), &
266 456 : out_d2cart=d2cart(:,:,:,ivol), out_eigvec=eigvec(:,:,:,ivol))
267 : end if
268 :
269 : call massmult_and_breaksym(gruns%cryst_vol(ivol)%natom, gruns%cryst_vol(ivol)%ntypat, &
270 912 : gruns%cryst_vol(ivol)%typat, gruns%ifc_vol(ivol)%amu, d2cart(:,:,:,ivol))
271 :
272 : !call zgemm('N','N',natom3,natom3,natom3,cone,d2cart(:,:,:,ivol),natom3,eigvec(:,:,:,ivol),natom3,czero,omat,natom3)
273 : !do nu=1,natom3
274 : ! write(std_out,*)"H|psi> - w**2 |psi>",maxval(abs(omat(:,:,nu) - wvols(nu,ivol) ** 2 * eigvec(:,:,nu,ivol)))
275 : !end do
276 : end do
277 :
278 : ! Compute dD(q)/dV with central finite difference.
279 57684 : dddv = zero
280 912 : do ivol=1,gruns%nvols
281 684 : fact = central_finite_diff(1, ivol, gruns%nvols)
282 115824 : if (fact /= zero) dddv = dddv + fact * d2cart(:,:,:,ivol)
283 : end do
284 57684 : dddv = dddv / gruns%delta_vol
285 :
286 : ! Compute -V0/(2w(q)**2) <u(q)|dD(q)/dq|u(q)>
287 228 : call zgemm('N','N',natom3,natom3,natom3,cone,dddv,natom3,eigvec(:,:,:,gruns%iv0),natom3,czero,omat,natom3)
288 2280 : do nu=1,natom3
289 2280 : if (abs(wvols(nu, gruns%iv0)) > tol12) then
290 2043 : dot = cg_zdotc(natom3, eigvec(1,1,nu, gruns%iv0), omat(1,1,nu))
291 : ! Must change sign if we have a purely imaginary solution i.e. z = iw
292 2043 : gvals(nu) = -sign(one, wvols(nu, gruns%iv0)) * gruns%v0 * dot(1) / (two * wvols(nu, gruns%iv0)**2)
293 : else
294 9 : gvals(nu) = zero
295 : end if
296 : end do
297 :
298 228 : end subroutine gruns_fourq
299 : !!***
300 :
301 : !----------------------------------------------------------------------
302 :
303 : !!****f* m_gruneisen/gruns_qpath
304 : !! NAME
305 : !! gruns_qpath
306 : !!
307 : !! FUNCTION
308 : !! Compute gruneisen parameters and group velocities on a q-path
309 : !! Write results to file.
310 : !!
311 : !! INPUTS
312 : !! prefix=Prefix for output files.
313 : !! qpath<kpath_t>=Object describing the q-path.
314 : !! ncid=netcdf file id.
315 : !! comm=MPI communicator.
316 : !!
317 : !! OUTPUT
318 : !! Only writing
319 : !!
320 : !! SOURCE
321 :
322 1 : subroutine gruns_qpath(gruns, prefix, qpath, ncid, comm)
323 :
324 : !Arguments ------------------------------------
325 : !scalars
326 : class(gruns_t),intent(in) :: gruns
327 : integer,intent(in) :: ncid,comm
328 : type(kpath_t),intent(in) :: qpath
329 : character(len=*),intent(in) :: prefix
330 :
331 : !Local variables-------------------------------
332 : !scalars
333 : integer,parameter :: master=0
334 : integer :: nprocs,my_rank,iqpt,ierr,ncerr,unt,iv0,nu,ii
335 : character(len=500) :: msg
336 : !arrays
337 1 : real(dp),allocatable :: gvals_qpath(:,:),wvols_qpath(:,:,:),dwdq_qpath(:,:,:)
338 1 : real(dp),allocatable :: phdispl_cart_qpath(:,:,:,:,:)
339 :
340 : ! ************************************************************************
341 :
342 1 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
343 :
344 81 : write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Gruneisen parameters along q-path ',ch10
345 1 : call wrtout(std_out, msg)
346 : !call wrtout(ab_out, msg)
347 :
348 1 : iv0 = gruns%iv0
349 2578 : ABI_CALLOC(wvols_qpath, (gruns%natom3, gruns%nvols, qpath%npts))
350 834 : ABI_CALLOC(gvals_qpath, (gruns%natom3, qpath%npts))
351 3075 : ABI_CALLOC(dwdq_qpath, (3, gruns%natom3, qpath%npts))
352 63086 : ABI_CALLOC(phdispl_cart_qpath, (2, gruns%natom3, gruns%natom3, gruns%nvols, qpath%npts))
353 :
354 84 : do iqpt=1,qpath%npts
355 83 : if (mod(iqpt, nprocs) /= my_rank) cycle ! mpi-parallelism
356 : call gruns_fourq(gruns, qpath%points(:,iqpt), wvols_qpath(:,:,iqpt), gvals_qpath(:,iqpt), &
357 84 : dwdq_qpath(:,:,iqpt), phdispl_cart_qpath(:,:,:,:,iqpt))
358 : end do
359 :
360 1 : call xmpi_sum(wvols_qpath, comm, ierr)
361 1 : call xmpi_sum(gvals_qpath, comm, ierr)
362 1 : call xmpi_sum(dwdq_qpath, comm, ierr)
363 1 : call xmpi_sum(phdispl_cart_qpath, comm, ierr)
364 :
365 : ! Write text files with phonon frequencies and gruneisen on the path.
366 1 : if (my_rank == master) then
367 1 : if (open_file(strcat(prefix, "_GRUNS_QPATH"), msg, newunit=unt, form="formatted", action="write") /= 0) then
368 1 : ABI_ERROR(msg)
369 : end if
370 1 : write(unt,'(a)')'# Phonon band structure, Gruneisen parameters and group velocity'
371 1 : write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
372 2 : call qpath%print([unt], pre="#")
373 : write(unt,'(5a)')&
374 1 : "# phfreq(mode=1) gruneisen(mode=1) velocity(mode=1) phfreq(mode=2) gruneisen(mode=2) velocity(mode=2) ..."
375 84 : do iqpt=1,qpath%npts
376 830 : do nu=1,gruns%natom3
377 : write(unt, "(3es17.8)", advance="no") &
378 3071 : wvols_qpath(nu, iv0, iqpt), gvals_qpath(nu, iqpt), sum(dwdq_qpath(1:3, nu, iqpt) ** 2)
379 : end do
380 84 : write(unt, "(a)")" "
381 : end do
382 1 : close(unt)
383 : end if
384 :
385 1 : if (my_rank == master .and. ncid /= nctk_noid) then
386 2 : ncerr = nctk_def_dims(ncid, [nctkdim_t("gruns_nqpath", qpath%npts)], defmode=.True.)
387 1 : NCF_CHECK(ncerr)
388 :
389 : ncerr = nctk_def_arrays(ncid, [ &
390 : ! q-points of the path
391 : nctkarr_t("gruns_qpath", "dp", "three, gruns_nqpath"), &
392 : ! gruneisen parameters on the path
393 : nctkarr_t("gruns_gvals_qpath", "dp", "number_of_phonon_modes, gruns_nqpath"), &
394 : ! phonon frequencies at the different volumes
395 : nctkarr_t("gruns_wvols_qpath", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqpath"), &
396 : ! group velocities at V0 in Cartesian coordinates.
397 : nctkarr_t("gruns_dwdq_qpath", "dp", "three, number_of_phonon_modes, gruns_nqpath"), &
398 : ! displacements for the different volumes.
399 : nctkarr_t("gruns_phdispl_cart_qpath", "dp", &
400 : "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqpath") &
401 6 : ])
402 1 : NCF_CHECK(ncerr)
403 :
404 : ! Write data.
405 1 : NCF_CHECK(nctk_set_datamode(ncid))
406 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qpath"), qpath%points))
407 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qpath"), gvals_qpath))
408 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qpath"), wvols_qpath))
409 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qpath"), dwdq_qpath))
410 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qpath"), phdispl_cart_qpath))
411 : end if
412 :
413 1 : ABI_FREE(wvols_qpath)
414 1 : ABI_FREE(gvals_qpath)
415 1 : ABI_FREE(dwdq_qpath)
416 1 : ABI_FREE(phdispl_cart_qpath)
417 :
418 1 : end subroutine gruns_qpath
419 : !!***
420 :
421 : !----------------------------------------------------------------------
422 :
423 : !!****f* m_gruneisen/gruns_qmesh
424 : !! NAME
425 : !! gruns_qmesh
426 : !!
427 : !! FUNCTION
428 : !! Compute gruneisen parameters and group velocities on a q-mesh.
429 : !! Save results to file.
430 : !!
431 : !! INPUTS
432 : !! prefix=Prefix for output files.
433 : !! dosdeltae=Step for the frequency mesh.
434 : !! ngqpt(3)=q-mesh divisions
435 : !! nshiftq=Number of shifts used to generated the ab-initio q-mesh.
436 : !! shiftq(3,nshiftq)=The shifts of the ab-initio q-mesh.
437 : !! ncid=netcdf file id.
438 : !! comm=MPI communicator
439 : !!
440 : !! OUTPUT
441 : !!
442 : !! SOURCE
443 :
444 1 : subroutine gruns_qmesh(gruns, prefix, dosdeltae, ngqpt, nshiftq, shiftq, ncid, comm)
445 :
446 : !Arguments ------------------------------------
447 : !scalars
448 : class(gruns_t),intent(in) :: gruns
449 : integer,intent(in) :: nshiftq,ncid,comm
450 : real(dp),intent(in) :: dosdeltae !,dossmear
451 : character(len=*),intent(in) :: prefix
452 : !arrays
453 : integer,intent(in) :: ngqpt(3)
454 : real(dp),intent(in) :: shiftq(3,nshiftq)
455 :
456 : !Local variables-------------------------------
457 : !scalars
458 : integer,parameter :: master=0,qptopt1=1,bcorr0=0
459 : integer :: nprocs,my_rank,iqibz,nqbz,nqibz,ierr,ii,nu,ncerr,nomega,cnt,unt,io
460 : real(dp) :: gavg,omega_min,omega_max,v2
461 1 : type(htetra_t) :: tetra
462 : character(len=500) :: msg
463 : !arrays
464 : integer :: qptrlatt(3,3)
465 1 : real(dp),allocatable :: gvals_qibz(:,:),wvols_qibz(:,:,:),dwdq_qibz(:,:,:)
466 1 : real(dp),allocatable :: qibz(:,:),qbz(:,:),wtq(:)
467 1 : real(dp),allocatable :: wdt(:,:),wdos(:,:),grdos(:,:),gr2dos(:,:),wibz(:),omega_mesh(:)
468 1 : real(dp),allocatable :: vdos(:,:),v2dos(:,:)
469 1 : real(dp),allocatable :: phdispl_cart_qibz(:,:,:,:,:)
470 :
471 : ! ************************************************************************
472 :
473 1 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
474 :
475 81 : write(msg,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of Gruneisen DOSes ',ch10
476 1 : call wrtout(std_out, msg)
477 :
478 : ! Generate the q-mesh by finding the IBZ and the corresponding weights.
479 4 : ABI_CHECK(all(ngqpt > 0), sjoin("invalid ngqpt:", ltoa(ngqpt)))
480 1 : qptrlatt = 0
481 4 : do ii=1,3
482 4 : qptrlatt(ii,ii) = ngqpt(ii)
483 : end do
484 :
485 : ! Get IBZ and BZ.
486 : call kpts_ibz_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptrlatt, qptopt1, nshiftq, shiftq, &
487 1 : nqibz, qibz, wtq, nqbz, qbz)
488 :
489 : ! Build tetrahedra
490 1 : tetra = tetra_from_kptrlatt(gruns%cryst_vol(gruns%iv0), qptopt1, qptrlatt, nshiftq, shiftq, nqibz, qibz, comm, msg, ierr)
491 1 : if (ierr /= 0) ABI_ERROR(msg)
492 :
493 4500 : ABI_CALLOC(wvols_qibz, (gruns%natom3, gruns%nvols, nqibz))
494 1454 : ABI_CALLOC(gvals_qibz, (gruns%natom3, nqibz))
495 5369 : ABI_CALLOC(dwdq_qibz, (3, gruns%natom3, nqibz))
496 110206 : ABI_CALLOC(phdispl_cart_qibz, (2, gruns%natom3, gruns%natom3, gruns%nvols, nqibz))
497 :
498 1 : gavg = zero
499 146 : do iqibz=1,nqibz
500 145 : if (mod(iqibz, nprocs) /= my_rank) cycle ! mpi-parallelism
501 : call gruns_fourq(gruns, qibz(:,iqibz), wvols_qibz(:,:,iqibz), gvals_qibz(:,iqibz), &
502 145 : dwdq_qibz(:,:,iqibz), phdispl_cart_qibz(:,:,:,:,iqibz))
503 1451 : gavg = gavg + wtq(iqibz) * sum(gvals_qibz(:,iqibz))
504 : end do
505 1 : gavg = gavg / gruns%natom3
506 :
507 1 : call xmpi_sum(gavg, comm, ierr)
508 1 : call xmpi_sum(wvols_qibz, comm, ierr)
509 1 : call xmpi_sum(gvals_qibz, comm, ierr)
510 1 : call xmpi_sum(dwdq_qibz, comm, ierr)
511 1 : call xmpi_sum(phdispl_cart_qibz, comm, ierr)
512 :
513 1 : omega_min = gruns%ifc_vol(gruns%iv0)%omega_minmax(1)
514 1 : omega_max = gruns%ifc_vol(gruns%iv0)%omega_minmax(2)
515 1 : nomega = nint((omega_max - omega_min) / dosdeltae) + 1
516 1 : nomega = max(6, nomega) ! Ensure Simpson integration will be ok
517 :
518 3 : ABI_MALLOC(omega_mesh, (nomega))
519 1 : omega_mesh = arth(omega_min, dosdeltae, nomega)
520 1 : omega_max = omega_mesh(nomega)
521 : !write(std_out,*)"hello",omega_min,omega_max,dosdeltae,(omega_max-omega_min) / (nomega-1)
522 3 : ABI_MALLOC(wibz, (nqibz))
523 3 : ABI_MALLOC(wdt, (nomega, 2))
524 3990 : ABI_CALLOC(wdos, (nomega, 2))
525 3990 : ABI_CALLOC(grdos, (nomega, 2))
526 3990 : ABI_CALLOC(gr2dos, (nomega, 2))
527 3990 : ABI_CALLOC(vdos, (nomega, 2))
528 3990 : ABI_CALLOC(v2dos, (nomega, 2))
529 :
530 : ! Compute DOSes.
531 1 : cnt = 0
532 146 : do iqibz=1,nqibz
533 1451 : do nu=1,gruns%natom3
534 1305 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi-parallelism
535 191835 : wibz = wvols_qibz(nu, gruns%iv0, :)
536 1305 : call tetra%get_onewk(iqibz,bcorr0,nomega,nqibz,wibz,omega_min,omega_max,one,wdt)
537 5205645 : wdt = wdt*wtq(iqibz)
538 5206950 : wdos = wdos + wdt
539 5206950 : grdos = grdos + wdt * gvals_qibz(nu,iqibz)
540 5206950 : gr2dos = gr2dos + wdt * gvals_qibz(nu,iqibz) ** 2
541 5220 : v2 = sum(dwdq_qibz(1:3,nu,iqibz) ** 2)
542 5206950 : vdos = vdos + wdt * sqrt(v2)
543 5207095 : v2dos = v2dos + wdt * v2
544 : end do
545 : end do
546 :
547 1 : call xmpi_sum(wdos, comm, ierr)
548 1 : call xmpi_sum(grdos, comm, ierr)
549 1 : call xmpi_sum(gr2dos, comm, ierr)
550 1 : call xmpi_sum(vdos, comm, ierr)
551 1 : call xmpi_sum(v2dos, comm, ierr)
552 :
553 1 : if (my_rank == master) then
554 1 : call wrtout(ab_out, sjoin(" Average Gruneisen parameter:", ftoa(gavg, fmt="f8.5")))
555 :
556 : ! Write text files with Gruneisen and DOSes.
557 1 : if (open_file(strcat(prefix, "_GRUNS_DOS"), msg, newunit=unt, form="formatted", action="write") /= 0) then
558 1 : ABI_ERROR(msg)
559 : end if
560 1 : write(unt,'(a)')'# Phonon density of states, Gruneisen DOS and phonon group velocity DOS'
561 1 : write(unt,'(a)')"# Energy in Hartree, DOS in states/Hartree"
562 1 : write(unt,'(a,i0)')'# Tetrahedron method with nqibz= ',nqibz
563 1 : write(unt,"(a,f8.5)")"# Average Gruneisen parameter:", gavg
564 : write(unt,'(5a)') &
565 1 : "# omega PH_DOS Gruns_DOS Gruns**2_DOS Vel_DOS Vel**2_DOS PH_IDOS Gruns_IDOS Gruns**2_IDOS Vel_IDOS Vel**2_IDOS"
566 1994 : do io=1,nomega
567 1993 : write(unt, "(11es17.8)")omega_mesh(io), &
568 1993 : wdos(io,1), grdos(io,1), gr2dos(io,1), vdos(io,1), v2dos(io,1), &
569 3987 : wdos(io,2), grdos(io,2), gr2dos(io,2), vdos(io,2), v2dos(io,2)
570 : end do
571 1 : close(unt)
572 : end if
573 :
574 : ! Write netcdf files.
575 1 : if (my_rank == master .and. ncid /= nctk_noid) then
576 : ncerr = nctk_def_dims(ncid, [ &
577 : nctkdim_t("gruns_nqibz", nqibz), nctkdim_t('gruns_nshiftq', nshiftq), &
578 4 : nctkdim_t('gruns_nomega', nomega)], defmode=.True.)
579 1 : NCF_CHECK(ncerr)
580 :
581 : ncerr = nctk_def_arrays(ncid, [ &
582 : ! q-point sampling in IBZ,
583 : nctkarr_t("gruns_qptrlatt", "int", "three, three"), &
584 : nctkarr_t("gruns_shiftq", "dp", "three, gruns_nshiftq"), &
585 : nctkarr_t("gruns_qibz", "dp", "three, gruns_nqibz"), &
586 : nctkarr_t("gruns_wtq", "dp", "gruns_nqibz"), &
587 : ! gruneisen parameters in IBZ
588 : ! phonon frequencies at the different volumes,
589 : ! group velocities at V0 in Cartesian coordinates.
590 : nctkarr_t("gruns_gvals_qibz", "dp", "number_of_phonon_modes, gruns_nqibz"), &
591 : nctkarr_t("gruns_wvols_qibz", "dp", "number_of_phonon_modes, gruns_nvols, gruns_nqibz"), &
592 : nctkarr_t("gruns_dwdq_qibz", "dp", "three, number_of_phonon_modes, gruns_nqibz"), &
593 : ! displacements for the different volumes.
594 : nctkarr_t("gruns_phdispl_cart_qibz", "dp", &
595 : "two, number_of_phonon_modes, number_of_phonon_modes, gruns_nvols, gruns_nqibz"), &
596 : ! DOSes and IDOSes
597 : nctkarr_t("gruns_omega_mesh", "dp", "gruns_nomega"), &
598 : nctkarr_t("gruns_wdos", "dp", "gruns_nomega, two"), &
599 : nctkarr_t("gruns_grdos", "dp", "gruns_nomega, two"), &
600 : nctkarr_t("gruns_gr2dos", "dp", "gruns_nomega, two"), &
601 : nctkarr_t("gruns_v2dos", "dp", "gruns_nomega, two"), &
602 : nctkarr_t("gruns_vdos", "dp", "gruns_nomega, two") &
603 15 : ])
604 1 : NCF_CHECK(ncerr)
605 :
606 : ! Write data.
607 1 : NCF_CHECK(nctk_set_datamode(ncid))
608 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qptrlatt"), qptrlatt))
609 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_shiftq"), shiftq))
610 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_qibz"), qibz))
611 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wtq"), wtq))
612 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gvals_qibz"), gvals_qibz))
613 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wvols_qibz"), wvols_qibz))
614 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_dwdq_qibz"), dwdq_qibz))
615 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_phdispl_cart_qibz"), phdispl_cart_qibz))
616 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_omega_mesh"), omega_mesh))
617 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_wdos"), wdos))
618 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_grdos"), grdos))
619 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_gr2dos"), gr2dos))
620 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_v2dos"), v2dos))
621 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_vdos"), vdos))
622 : end if
623 :
624 1 : ABI_FREE(qibz)
625 1 : ABI_FREE(wtq)
626 1 : ABI_FREE(qbz)
627 1 : ABI_FREE(wvols_qibz)
628 1 : ABI_FREE(gvals_qibz)
629 1 : ABI_FREE(dwdq_qibz)
630 1 : ABI_FREE(phdispl_cart_qibz)
631 1 : ABI_FREE(omega_mesh)
632 1 : ABI_FREE(wibz)
633 1 : ABI_FREE(wdt)
634 1 : ABI_FREE(wdos)
635 1 : ABI_FREE(grdos)
636 1 : ABI_FREE(gr2dos)
637 1 : ABI_FREE(v2dos)
638 1 : ABI_FREE(vdos)
639 :
640 1 : call tetra%free()
641 :
642 3 : end subroutine gruns_qmesh
643 : !!***
644 :
645 : !----------------------------------------------------------------------
646 :
647 : !!****f* m_gruneisen/gruns_free
648 : !! NAME
649 : !! gruns_free
650 : !!
651 : !! FUNCTION
652 : !! Free dynamic memory.
653 : !!
654 : !! SOURCE
655 :
656 1 : subroutine gruns_free(gruns)
657 :
658 : !Arguments ------------------------------------
659 : class(gruns_t),intent(inout) :: gruns
660 :
661 : !Local variables-------------------------------
662 : integer :: ii
663 : ! ************************************************************************
664 :
665 1 : if (allocated(gruns%ifc_vol)) then
666 4 : do ii=1,size(gruns%cryst_vol)
667 4 : call gruns%cryst_vol(ii)%free()
668 : end do
669 4 : ABI_FREE(gruns%cryst_vol)
670 : end if
671 :
672 1 : if (allocated(gruns%ddb_vol)) then
673 4 : do ii=1,size(gruns%ddb_vol)
674 4 : call gruns%ddb_vol(ii)%free()
675 : end do
676 4 : ABI_FREE(gruns%ddb_vol)
677 : end if
678 :
679 1 : if (allocated(gruns%ifc_vol)) then
680 4 : do ii=1,size(gruns%ifc_vol)
681 4 : call gruns%ifc_vol(ii)%free()
682 : end do
683 4 : ABI_FREE(gruns%ifc_vol)
684 : end if
685 :
686 1 : end subroutine gruns_free
687 : !!***
688 :
689 : !----------------------------------------------------------------------
690 :
691 : !!****f* m_gruneisen/gruns_anaddb
692 : !! NAME
693 : !! gruns_anaddb
694 : !!
695 : !! FUNCTION
696 : !! Driver routine called in anaddb to compute Gruneisen parameters.
697 : !!
698 : !! INPUTS
699 : !! inp<anaddb_dataset_type: :: anaddb input variables.
700 : !! prefix=Prefix for output files
701 : !! comm=MPI communicator
702 : !!
703 : !! OUTPUT
704 : !! Only writing.
705 : !!
706 : !! SOURCE
707 :
708 1 : subroutine gruns_anaddb(inp, comm)
709 :
710 : !Arguments ------------------------------------
711 : integer,intent(in) :: comm
712 : type(anaddb_dataset_type),intent(inout) :: inp
713 :
714 : !Local variables-------------------------------
715 : !scalars
716 : integer,parameter :: master=0
717 : integer :: ii,nprocs,my_rank,ncid,iv0, ncerr
718 : real(dp) :: cpu,wall,gflops
719 : type(gruns_t),target :: gruns
720 1 : type(kpath_t) :: qpath
721 : ! ************************************************************************
722 :
723 1 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
724 :
725 1 : ABI_CHECK(inp%ifcflag == 1, "Gruneisen requires ifcflag == 1")
726 :
727 1 : call cwtime(cpu, wall, gflops, "start")
728 :
729 1 : gruns = gruns_new(inp%gruns_ddbs, inp, comm)
730 1 : iv0 = gruns%iv0
731 :
732 1 : ncid = nctk_noid
733 1 : if (my_rank == master) then
734 1 : NCF_CHECK_MSG(nctk_open_create(ncid, strcat(inp%filename_output, "_GRUNS.nc"), xmpi_comm_self), "Creating _GRUNS.nc")
735 :
736 : ! Write structure corresponding to iv0
737 1 : NCF_CHECK(gruns%cryst_vol(iv0)%ncwrite(ncid))
738 :
739 : ! Add important dimensions and additional metadata.
740 : ncerr = nctk_def_dims(ncid, [ &
741 : nctkdim_t("gruns_nvols", gruns%nvols), &
742 3 : nctkdim_t('number_of_phonon_modes', gruns%natom3)], defmode=.True.)
743 1 : NCF_CHECK(ncerr)
744 2 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"])
745 1 : NCF_CHECK(ncerr)
746 :
747 : ! Add lattice parameters and positions of the `gruns_nvols` structures so
748 : ! that we can easily reconstruct structure objects in AbiPy.
749 : ncerr = nctk_def_arrays(ncid, [ &
750 : nctkarr_t("gruns_rprimd", "dp", "three, three, gruns_nvols"), &
751 : nctkarr_t("gruns_xred", "dp", "three, number_of_atoms, gruns_nvols") &
752 3 : ])
753 1 : NCF_CHECK(ncerr)
754 :
755 3 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "gruns_iv0"], [iv0], datamode=.True.)
756 1 : NCF_CHECK(ncerr)
757 :
758 4 : do ii=1,gruns%nvols
759 12 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_rprimd"), gruns%cryst_vol(ii)%rprimd, start=[1,1,ii]))
760 13 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gruns_xred"), gruns%cryst_vol(ii)%xred, start=[1,1,ii]))
761 : end do
762 :
763 : !call phonons_ncwrite(ncid,natom,nfineqpath,save_qpoints,weights,save_phfrq,save_phdispl_cart)
764 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'atomic_mass_units'), ddb%amu))
765 : end if
766 :
767 : ! Compute gruneisen parameters on the q-mesh.
768 4 : if (all(inp%ng2qpt /= 0)) then
769 1 : call gruns_qmesh(gruns, inp%filename_output, inp%dosdeltae, inp%ng2qpt, 1, inp%q2shft, ncid, comm)
770 : else
771 0 : ABI_WARNING("Cannot compute Gruneisen parameters on q-mesh because ng2qpt == 0")
772 : end if
773 :
774 : ! Compute gruneisen on the q-path.
775 1 : if (inp%nqpath /= 0) then
776 1 : call qpath%init(inp%qpath, gruns%cryst_vol(iv0)%gprimd, inp%ndivsm)
777 1 : call gruns_qpath(gruns, inp%filename_output, qpath, ncid, comm)
778 1 : call qpath%free()
779 : else
780 0 : ABI_WARNING("Cannot compute Gruneisen parameters on q-path because nqpath == 0")
781 : end if
782 :
783 : ! Compute speed of sound for V0.
784 1 : if (inp%vs_qrad_tolkms(1) > zero) then
785 1 : call gruns%ifc_vol(iv0)%speedofsound(gruns%cryst_vol(iv0), inp%vs_qrad_tolkms, ncid, comm)
786 : end if
787 :
788 : ! Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities)
789 1 : if (my_rank == master .and. inp%nph2l /= 0 .and. inp%ifcflag == 1) then
790 0 : call gruns%ifc_vol(iv0)%calcnwrite_nana_terms(gruns%cryst_vol(iv0), inp%nph2l, inp%qph2l, inp%qnrml2, ncid)
791 : end if
792 :
793 : if (my_rank == master) then
794 1 : NCF_CHECK(nf90_close(ncid))
795 : end if
796 :
797 1 : call gruns_free(gruns)
798 1 : call cwtime_report("gruns_anaddb", cpu, wall, gflops)
799 :
800 1 : end subroutine gruns_anaddb
801 : !!***
802 :
803 : !----------------------------------------------------------------------
804 :
805 0 : end module m_gruneisen
|
