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
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
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
|
