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