Branch data Line data Source code
1 : : !!****m*ABINIT/m_anaddb_driver
2 : : !! NAME
3 : : !! m_anaddb_driver
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2014-2026 ABINIT group (GA)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_anaddb_driver
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_xmpi
26 : : use m_errors
27 : : use m_nctk
28 : : use netcdf
29 : : use m_ifc
30 : : use m_ddb_hdr
31 : : use m_phonons
32 : : use m_supercell
33 : : use m_raman
34 : :
35 : : use m_fstrings, only : strcat
36 : : use m_crystal, only : crystal_t
37 : : use m_kpts, only : qmesh_from_brav_ngqpt
38 : : use m_anaddb_dataset, only : anaddb_dataset_type
39 : : use m_ddb, only : ddb_type, asrq0_t
40 : : use m_dynmat, only : gtdyn9, dfpt_phfrq, dfpt_prtph
41 : : use m_ddb_interpolate, only : ddb_interpolate
42 : : use m_harmonic_thermo, only : harmonic_thermo
43 : : use m_elphon, only : elphon
44 : : use m_thmeig, only : thmeig
45 : : use m_relaxpol, only : relaxpol
46 : : use m_ddb_diel, only : ddb_diel
47 : : use m_ddb_elast, only : ddb_elast
48 : : use m_ddb_piezo, only : ddb_piezo
49 : : use m_ddb_internalstr, only : ddb_internalstr
50 : : use m_gruneisen, only : gruns_anaddb
51 : : use m_ddb_flexo, only : ddb_flexo, ddb_phi1
52 : : use m_lwf, only : run_lattice_wannier
53 : :
54 : : implicit none
55 : :
56 : : private
57 : :
58 : : public:: anaddb_driver_type
59 : :
60 : : !----------------------------------------------------------------------
61 : :
62 : : !!****t* m_anaddb_driver/anaddb_driver_type
63 : : !! NAME
64 : : !! anaddb_driver_type
65 : : !!
66 : : !! FUNCTION
67 : : !! The anaddb_driver_type structured datatype
68 : : !! Main subroutines run by anaddb.
69 : : !!
70 : : !! SOURCE
71 : :
72 : : type anaddb_driver_type
73 : :
74 : : integer:: natom
75 : : integer:: msize
76 : : integer:: mpert
77 : :
78 : : logical:: do_ifc=.false.
79 : : logical:: do_electric_tensors=.false.
80 : : logical:: do_dielectric_q0=.false.
81 : : logical:: do_dielectric_nonana=.false.
82 : : logical:: do_phonon_dos=.false.
83 : : logical:: do_phonon_bs=.false.
84 : :
85 : : real(dp):: epsinf(3, 3)
86 : : real(dp):: dchide(3,3,3)
87 : : real(dp):: dielt_rlx(3, 3)
88 : : real(dp):: elast(6, 6)
89 : :
90 : : real(dp), allocatable:: d2cart(:,:)
91 : : ! d2cart(2,msize)
92 : :
93 : : real(dp), allocatable:: displ(:)
94 : : ! displ(2*3*natom*3*natom)
95 : :
96 : : real(dp), allocatable:: phfrq(:)
97 : : ! phfrq(3*natom)
98 : :
99 : : real(dp), allocatable:: instrain(:,:)
100 : : ! instrain(3*natom,6)
101 : :
102 : : real(dp), allocatable:: dchidt(:,:,:,:)
103 : : ! dchidt(natom,3,3,3)
104 : :
105 : : real(dp), allocatable:: fact_oscstr(:,:,:)
106 : : ! fact_oscstr(2,3,3*natom)
107 : :
108 : : real(dp), allocatable:: zeff(:,:,:)
109 : : ! zeff(3,3,natom)
110 : :
111 : : real(dp), allocatable:: qdrp_cart(:,:,:,:)
112 : : ! qdrp_cart(3,3,3,natom)
113 : :
114 : : real(dp), allocatable:: dcdq(:,:,:,:,:)
115 : : ! dcdq(3,natom,3,natom,3)
116 : :
117 : : real(dp), allocatable:: dcdqdq(:,:,:,:,:)
118 : : ! dcdqdq(3,natom,3,3,3)
119 : :
120 : : contains
121 : :
122 : : procedure :: init => anaddb_driver_init
123 : : ! Initialize object
124 : :
125 : : procedure :: free => anaddb_driver_free
126 : : ! Free memory
127 : :
128 : : procedure :: open_write_nc => anaddb_driver_open_write_nc
129 : : ! Open netcdf output and write some info
130 : :
131 : : procedure :: electric_tensors => anaddb_driver_electric_tensors
132 : : ! Compute dielectric tensor, born effective charges, and quadrupoles
133 : :
134 : : procedure :: structural_response => anaddb_driver_structural_response
135 : : ! Structural response at fixed polarization
136 : :
137 : : procedure :: susceptibilities => anaddb_driver_susceptibilities
138 : : ! Compute non-linear optical susceptibilities and first-order derivatives
139 : :
140 : : procedure :: interatomic_force_constants => anaddb_driver_interatomic_force_constants
141 : : ! Compute the interatomic force constants from a ddb
142 : :
143 : : procedure :: phdos => anaddb_driver_phdos
144 : : ! Compute phonon density of states
145 : :
146 : : procedure :: thermal_supercell => anaddb_driver_thermal_supercell
147 : : ! Thermal supercell calculation
148 : :
149 : : procedure :: harmonic_thermo => anaddb_driver_harmonic_thermo
150 : : ! Phonon density of states and thermodynamical properties
151 : :
152 : : procedure :: dielectric_q0 => anaddb_driver_dielectric_q0
153 : : ! Dielectric tensor and related properties
154 : :
155 : : procedure :: nonlinear_response => anaddb_driver_nonlinear_response
156 : : ! Non-linear response: electrooptic and Raman
157 : :
158 : : procedure :: dielectric_nonana => anaddb_driver_dielectric_nonana
159 : : ! Non-analyticity in the dielectric matrix and raman susceptibility
160 : :
161 : : procedure :: internal_strain => anaddb_driver_internal_strain
162 : : ! Internal strain
163 : :
164 : : procedure :: elastic_tensor => anaddb_driver_elastic_tensor
165 : : ! Elastic tensor
166 : :
167 : : procedure :: piezoelectric_tensor => anaddb_driver_piezoelectric_tensor
168 : : ! Piezoelectric tensor
169 : :
170 : : procedure :: flexoelectric_tensor => anaddb_driver_flexoelectric_tensor
171 : : ! Flexoelectric tensor
172 : :
173 : : procedure :: get_dcdq => anaddb_driver_get_dcdq
174 : : ! Get IFCs derivatives for the imposition of rotational invariance
175 : :
176 : : procedure :: lattice_wannier => anaddb_driver_lattice_wannier
177 : : ! Construct the Lattice Wannier functions
178 : :
179 : : procedure :: convertdim_dielt => anaddb_driver_convertdim_dielt
180 : : ! Extract low-dimmensional dielectric response from periodic calculations
181 : :
182 : : end type anaddb_driver_type
183 : : !!***
184 : :
185 : : contains
186 : :
187 : : !!****f* m_anaddb_driver/anaddb_driver_init
188 : : !! NAME
189 : : !! anaddb_driver_init
190 : : !!
191 : : !! FUNCTION
192 : : !! Initialize object
193 : : !!
194 : : !! INPUTS
195 : : !!
196 : : !! OUTPUT
197 : : !!
198 : : !! SOURCE
199 : :
200 : 81 : subroutine anaddb_driver_init(driver, dtset)
201 : :
202 : : !Arguments -------------------------------
203 : : class(anaddb_driver_type), intent(inout):: driver
204 : : type(anaddb_dataset_type), intent(in):: dtset
205 : : ! ************************************************************************
206 : :
207 : : ! Set control flags
208 [ + + ]: 81 : if (dtset%ifcflag == 1) then
209 : 55 : driver%do_ifc = .true.
210 : : end if
211 : :
212 : : if (dtset%ifcflag /= 0 .or. dtset%dieflag /= 0 &
213 : : & .or. dtset%dipdip /= 0 .or. dtset%dipquad /= 0 .or. dtset%quadquad /= 0 &
214 : : & .or. dtset%nph2l /= 0 .or. dtset%nlflag /= 0 &
215 : : & .or. dtset%piezoflag /= 0 .or. dtset%flexoflag /= 0 &
216 [ + + + + : 81 : & .or. dtset%polflag /= 0) then
+ - + - +
- + - + -
+ - - + ]
217 : 77 : driver%do_electric_tensors = .true.
218 : : end if
219 : :
220 : : if ((dtset%dieflag /= 0 .and. dtset%dieflag /= 2) &
221 : : & .or. dtset%nph2l /= 0 .or. dtset%nlflag == 1 &
222 [ + + + + : 81 : & .or. dtset%piezoflag /= 0) then
+ - + + ]
223 : 27 : driver%do_dielectric_q0 = .true.
224 : : end if
225 : :
226 [ + + ]: 81 : if (dtset%nph2l /= 0) then
227 : 24 : driver%do_dielectric_nonana = .true.
228 : : end if
229 : :
230 [ + + + + ]: 81 : if (dtset%ifcflag == 1 .and. any(dtset%prtdos==[1, 2])) then
231 : 11 : driver%do_phonon_dos = .true.
232 : : end if
233 : :
234 [ + + + + ]: 81 : if (dtset%nph1l /= 0 .or. dtset%nqpath /= 0) then
235 : 66 : driver%do_phonon_bs = .true.
236 : : end if
237 : :
238 [ + + ]: 81 : if (dtset%gruns_nddbs /= 0) then
239 : 1 : driver%do_ifc = .false.
240 : 1 : driver%do_electric_tensors = .false.
241 : 1 : driver%do_dielectric_q0 = .false.
242 : 1 : driver%do_dielectric_nonana = .false.
243 : 1 : driver%do_phonon_bs = .false.
244 : 1 : driver%do_phonon_dos = .false.
245 : : end if
246 : :
247 : : ! Copy dimensions
248 : 81 : driver%natom = dtset%natom
249 : 81 : driver%msize = dtset%msize
250 : 81 : driver%mpert = dtset%mpert
251 : :
252 : : ! Allocate memory
253 [ - + + - : 243 : ABI_MALLOC(driver%d2cart, (2, driver%msize))
+ - + - -
+ - + -
+ ]
254 [ + - - + : 243 : ABI_MALLOC(driver%displ, (2*3*driver%natom*3*driver%natom))
+ - + - +
- - + - +
- + ]
255 [ + - - + : 243 : ABI_MALLOC(driver%phfrq, (3*driver%natom))
+ - + - +
- - + - +
- + ]
256 [ - + + - : 324 : ABI_MALLOC(driver%instrain, (3*driver%natom, 6))
+ - + - +
- - + - +
- + ]
257 [ - + + - : 486 : ABI_MALLOC(driver%dcdq, (3, driver%natom, 3, driver%natom, 3))
+ - - + +
- + - + -
+ - - + -
+ - + ]
258 [ - + + - : 405 : ABI_MALLOC(driver%dcdqdq, (3, driver%natom, 3, 3, 3))
+ - + - +
- + - - +
- + - + ]
259 [ + + + + : 50820 : driver%dcdq = zero
+ + + + +
+ ]
260 [ + + + + : 31320 : driver%dcdqdq = zero
+ + + + +
+ ]
261 : : ! Electric tensors
262 [ - + + - : 243 : ABI_MALLOC(driver%zeff, (3, 3, driver%natom))
+ - + - -
+ - + -
+ ]
263 [ - + + - : 243 : ABI_MALLOC(driver%qdrp_cart, (3, 3, 3, driver%natom))
+ - + - -
+ - + -
+ ]
264 : :
265 : : ! oscillator strength and Lyddane-Sachs-Teller relation
266 [ - + + - : 162 : ABI_MALLOC(driver%fact_oscstr, (2, 3, 3*driver%natom))
- + - + -
+ ]
267 : :
268 : : ! Susceptibilities
269 [ + + ]: 81 : if (dtset%nlflag > 0) then
270 [ - + + - : 8 : ABI_MALLOC(driver%dchidt, (driver%natom, 3, 3, 3))
- + - + -
+ ]
271 : : end if
272 : :
273 : 81 : end subroutine anaddb_driver_init
274 : : !!***
275 : :
276 : :
277 : : !!****f* m_anaddb_driver/anaddb_driver_free
278 : : !! NAME
279 : : !! anaddb_driver_free
280 : : !!
281 : : !! FUNCTION
282 : : !! Deallocate memory
283 : : !!
284 : : !! INPUTS
285 : : !!
286 : : !! OUTPUT
287 : : !!
288 : : !! SOURCE
289 : :
290 : 81 : subroutine anaddb_driver_free(driver)
291 : :
292 : : !Arguments -------------------------------
293 : : class(anaddb_driver_type), intent(inout):: driver
294 : : ! ************************************************************************
295 : :
296 [ + - ]: 81 : ABI_SFREE(driver%zeff)
297 [ + - ]: 81 : ABI_SFREE(driver%qdrp_cart)
298 [ + - ]: 81 : ABI_SFREE(driver%d2cart)
299 [ + - ]: 81 : ABI_SFREE(driver%displ)
300 [ + - ]: 81 : ABI_SFREE(driver%phfrq)
301 [ + + ]: 81 : ABI_SFREE(driver%dchidt)
302 [ + - ]: 81 : ABI_SFREE(driver%instrain)
303 [ + - ]: 81 : ABI_SFREE(driver%fact_oscstr)
304 [ + - ]: 81 : ABI_SFREE(driver%dcdq)
305 [ + - ]: 81 : ABI_SFREE(driver%dcdqdq)
306 : :
307 : 81 : end subroutine anaddb_driver_free
308 : : !!***
309 : :
310 : : !!****f* m_anaddb_driver/anaddb_driver_open_write_nc
311 : : !! NAME
312 : : !! anaddb_driver_open_write_nc
313 : : !!
314 : : !! FUNCTION
315 : : !! Open anaddb netcdf output and define dimensions.
316 : : !!
317 : : !! INPUTS
318 : : !!
319 : : !! OUTPUT
320 : : !!
321 : : !! SOURCE
322 : :
323 : 81 : subroutine anaddb_driver_open_write_nc(driver, ana_ncid, dtset, crystal, comm)
324 : :
325 : : !Arguments -------------------------------
326 : : class(anaddb_driver_type), intent(inout):: driver
327 : : type(anaddb_dataset_type), intent(in):: dtset
328 : : type(crystal_t),intent(in):: crystal
329 : : integer,intent(in):: comm
330 : : integer, intent(out):: ana_ncid
331 : :
332 : : !Local variables -------------------------------
333 : : integer, parameter:: master = 0
334 : : integer:: natom,lenstr
335 : : integer:: ncerr
336 : : integer:: my_rank
337 : : ! ************************************************************************
338 : :
339 : 81 : my_rank = xmpi_comm_rank(comm)
340 : :
341 : 81 : natom = driver%natom
342 : 81 : lenstr = dtset%lenstr
343 : :
344 : : ! Open the netcdf file that will contain the anaddb results
345 : 81 : ana_ncid = nctk_noid
346 [ + - ]: 81 : if (my_rank == master) then
347 [ - + - - ]: 81 : NCF_CHECK_MSG(nctk_open_create(ana_ncid, trim(dtset%prefix_outdata)//"_anaddb.nc", xmpi_comm_self), "Creating anaddb.nc")
348 : : ncerr = nctk_def_dims(ana_ncid, [ &
349 : : nctkdim_t('number_of_atoms', natom), &
350 : : nctkdim_t('natom3', 3*natom), &
351 : : nctkdim_t('number_of_phonon_modes', 3*natom), &
352 : : nctkdim_t('anaddb_input_len', lenstr) &
353 [ + + ]: 405 : ], defmode=.True.)
354 [ - + ]: 81 : NCF_CHECK(ncerr)
355 : : ncerr = nctk_def_arrays(ana_ncid, [ &
356 : : nctkarr_t("anaddb_input_string", "char", "anaddb_input_len") &
357 [ + + ]: 162 : ])
358 [ - + ]: 81 : NCF_CHECK(ncerr)
359 : : !NCF_CHECK(nctk_defnwrite_ivars(ana_ncid, ["anaddb_version"], [1]))
360 : :
361 : : ncerr = nctk_def_iscalars(ana_ncid, [character(len = nctk_slen) :: &
362 [ + + ]: 567 : "asr", "chneut", "dipdip", "symdynmat", "dipquad", "quadquad"])
363 [ - + ]: 81 : NCF_CHECK(ncerr)
364 : :
365 [ - + ]: 81 : NCF_CHECK(nctk_set_datamode(ana_ncid))
366 : 81 : ncerr = nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "anaddb_input_string"), dtset%input_string(:lenstr))
367 [ - + ]: 81 : NCF_CHECK(ncerr)
368 [ - + ]: 81 : NCF_CHECK(crystal%ncwrite(ana_ncid))
369 : :
370 : : ncerr = nctk_write_iscalars(ana_ncid, [character(len = nctk_slen) :: &
371 : : "asr", "chneut", "dipdip", "symdynmat", "dipquad", "quadquad"], &
372 : : [dtset%asr, dtset%chneut, dtset%dipdip, &
373 [ + + + + ]: 1053 : dtset%symdynmat, dtset%dipquad, dtset%quadquad])
374 [ - + ]: 81 : NCF_CHECK(ncerr)
375 : : end if
376 : :
377 : 81 : end subroutine anaddb_driver_open_write_nc
378 : : !!***
379 : :
380 : : !!****f* m_anaddb_driver/anaddb_driver_electric_tensors
381 : : !! NAME
382 : : !! anaddb_driver_electric_tensors
383 : : !!
384 : : !! FUNCTION
385 : : !! get Dielectric tensor, born effective charges, and quadrupole tensor,
386 : : !! and write them to netcdf output.
387 : : !!
388 : : !! INPUTS
389 : : !!
390 : : !! OUTPUT
391 : : !!
392 : : !! SOURCE
393 : :
394 : 76 : subroutine anaddb_driver_electric_tensors(driver, dtset, crystal, ddb, ddb_lw, ddb_hdr, ana_ncid, comm)
395 : :
396 : : !Arguments -------------------------------
397 : : class(anaddb_driver_type), intent(inout):: driver
398 : : type(anaddb_dataset_type), intent(inout):: dtset
399 : : type(crystal_t), intent(in):: crystal
400 : : type(ddb_type), intent(inout):: ddb
401 : : type(ddb_type), intent(inout):: ddb_lw
402 : : type(ddb_hdr_type), intent(in):: ddb_hdr
403 : : integer, intent(in):: ana_ncid
404 : : integer, intent(in):: comm
405 : :
406 : : !Local variables -------------------------------
407 : : integer, parameter:: master = 0
408 : : integer:: my_rank
409 : : integer:: ii
410 : : integer:: iblok, iblok_quadrupoles, iblok_epsinf
411 : : integer:: ncerr
412 : : integer:: lwsym
413 : : character(len = 500):: msg
414 : : integer:: units(2)
415 : : ! ************************************************************************
416 : :
417 [ + + ]: 228 : units = [std_out, ab_out]
418 : :
419 : : ! Get Quadrupole tensor
420 : 76 : iblok_quadrupoles = 0
421 [ + + + + : 9596 : driver%qdrp_cart = zero
+ + + + ]
422 [ + + ]: 76 : if (ddb_hdr%has_d3E_lw) then
423 [ + - + + ]: 648 : write(msg, '(2a, (80a), 2a)') ch10, ('=',ii = 1, 80)
424 : 8 : call wrtout(units, msg)
425 : 8 : lwsym = 1
426 : 8 : iblok_quadrupoles = ddb_lw%get_quadrupoles(ddb_hdr%ddb_version, lwsym, BLKTYP_d3E_lw, driver%qdrp_cart)
427 : : end if
428 : :
429 : : ! Get the electronic dielectric tensor (epsinf) and Born effective charges (zeff)
430 : : ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
431 : 76 : iblok = ddb%get_dielt_zeff(crystal, dtset%rfmeth, dtset%chneut, dtset%selectz, driver%epsinf, driver%zeff)
432 : :
433 : : ! Try to get epsinf, in case just the DDE are present
434 [ + + ]: 76 : if (iblok == 0) then
435 : 32 : iblok_epsinf = ddb%get_dielt(dtset%rfmeth, driver%epsinf)
436 : : end if
437 : :
438 : : !if (iblok_epsinf == 0) then
439 : : !GA: Not the greatest way of checking
440 [ + + + - : 76 : if (driver%epsinf(1, 1)==one .and. driver%epsinf(2, 2)==one .and. driver%epsinf(3, 3)==one) then
+ - ]
441 : 33 : write(msg, '(5a)') ch10, &
442 : 33 : ' The DDB file does not contain the derivatives w.r.t. electric field perturbation. ',ch10, &
443 : 66 : ' The program will continue by setting the electronic dielectric tensor to 1. ',ch10
444 : : ! call wrtout([ab_out], msg)
445 : 33 : ABI_WARNING(msg)
446 : : end if
447 : :
448 : : !**********************************************************************
449 : : ! Write Dielectric tensor, born effective charges, and quadrupoles to netcdf output.
450 : :
451 : 76 : my_rank = xmpi_comm_rank(comm)
452 [ + - ]: 76 : if (my_rank == master) then
453 : : ncerr = nctk_def_arrays(ana_ncid, [&
454 : : nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'), &
455 : : nctkarr_t('quadrupoles_cart', "dp", 'three, three, three, number_of_atoms'), &
456 : : nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms")],&
457 [ + + ]: 304 : defmode=.True.)
458 [ - + ]: 76 : NCF_CHECK(ncerr)
459 : :
460 [ - + ]: 76 : NCF_CHECK(nctk_set_datamode(ana_ncid))
461 [ - + ]: 76 : NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, 'emacro_cart'), driver%epsinf))
462 [ - + ]: 76 : NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, 'quadrupoles_cart'), driver%qdrp_cart))
463 [ - + ]: 76 : NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, 'becs_cart'), driver%zeff))
464 : : end if
465 : :
466 : 76 : end subroutine anaddb_driver_electric_tensors
467 : : !!***
468 : :
469 : : !!****f* m_anaddb_driver/anaddb_driver_structural_response
470 : : !! NAME
471 : : !! anaddb_driver_structural_response
472 : : !!
473 : : !! FUNCTION
474 : : !! Compute structural response at fixed polarization
475 : : !!
476 : : !! INPUTS
477 : : !!
478 : : !! OUTPUT
479 : : !!
480 : : !! SOURCE
481 : :
482 : 3 : subroutine anaddb_driver_structural_response(driver, dtset, crystal, ddb)
483 : :
484 : : !Arguments -------------------------------
485 : : class(anaddb_driver_type), intent(inout):: driver
486 : : type(anaddb_dataset_type), intent(in):: dtset
487 : : type(crystal_t), intent(in):: crystal
488 : : type(ddb_type), intent(inout):: ddb
489 : :
490 : : !Local variables -------------------------------
491 : : integer:: iblok
492 : : integer:: msize
493 : : real(dp):: etotal
494 : : character(len = 500):: msg
495 : : integer:: rfelfd(4), rfphon(4), rfstrs(4)
496 : : real(dp):: red_ptot(3)
497 : : real(dp):: pel(3)
498 : : real(dp):: strten(6)
499 : : real(dp) :: targetpol(3)
500 : : real(dp):: qphnrm(3), qphon(3, 3)
501 : 3 : integer, allocatable:: d2flg(:)
502 : 3 : real(dp), allocatable:: gred(:,:)
503 : : ! ************************************************************************
504 : :
505 : 3 : msize = dtset%msize
506 [ - + + - : 9 : ABI_MALLOC(d2flg, (msize))
+ - + - -
+ - + ]
507 : :
508 : : ! Look for the Gamma Block in the DDB
509 [ + + ]: 12 : qphon(:,1)=zero
510 : 3 : qphnrm(1)=zero
511 [ + + ]: 9 : rfphon(1:2)=1
512 [ + + ]: 9 : rfelfd(1:2)=2
513 [ + + ]: 9 : rfstrs(1:2)=0
514 : :
515 : : !write(std_out,*)"ddb%mpert",ddb%mpert
516 : 3 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, dtset%rfmeth)
517 : : !iblok = ddb%get_dielt_zeff(crystal, dtset%rfmeth, dtset%chneut, dtset%selectz, driver%epsinf, driver%zeff)
518 : :
519 [ + + ]: 3 : if(iblok /= 0)then
520 : : ! Save the second-order derivatives
521 [ + + + + ]: 17984 : driver%d2cart(1:2, 1:msize) = ddb%val(1:2, 1:msize, iblok)
522 [ + + ]: 5996 : d2flg(1:msize) = ddb%flg(1:msize, iblok)
523 : :
524 : : else
525 : : ! the gamma blok has not been found
526 [ + - + - ]: 1 : if (dtset%relaxat == 0 .and. dtset%relaxstr == 0) then
527 : : ! The gamma blok is not needed
528 [ + + + + ]: 6076 : driver%d2cart(1:2, 1:msize)=zero
529 [ + + ]: 2026 : d2flg(1:msize)=1
530 : : else
531 : : ! There is a problem !
532 : : write(msg, '(7a)' )&
533 : 0 : 'The dynamical matrix at Gamma is needed, in order to perform ',ch10, &
534 : 0 : "relaxation at constant polarisation (Na Sai's method)",ch10, &
535 : 0 : 'However, this was not found in the DDB.',ch10, &
536 : 0 : 'Action: complete your DDB with the dynamical matrix at Gamma.'
537 : 0 : ABI_ERROR(msg)
538 : : end if
539 : : end if ! iblok not found
540 : :
541 : : ! Extract the block with the total energy
542 [ - + ]: 3 : if (ddb%get_etotal(etotal) == 0) then
543 : 0 : ABI_ERROR("DDB file does not contain GS etotal")
544 : : end if
545 : :
546 : : ! Extract the polarizability
547 : 3 : iblok = ddb%get_pel(pel, dtset%relaxat, dtset%relaxstr)
548 : :
549 : : ! Extract the forces
550 [ - + ]: 3 : iblok = ddb%get_gred(gred, dtset%relaxat, dtset%relaxstr)
551 : :
552 : : ! Extract the stress tensor
553 : 3 : iblok = ddb%get_strten(strten, dtset%relaxat, dtset%relaxstr)
554 : :
555 : : ! when called from here red_ptot is not set ! So set it to zero
556 : 3 : red_ptot(:) = zero
557 : :
558 [ + + ]: 12 : targetpol(:) = dtset%targetpol
559 : :
560 : : call relaxpol(crystal, d2flg, driver%d2cart, etotal, gred, dtset%iatfix, &
561 : : & ab_out, dtset%istrfix, dtset%mpert, dtset%msize, dtset%natfix, crystal%natom, &
562 : : & dtset%nstrfix, pel, red_ptot, dtset%relaxat, dtset%relaxstr, &
563 : 3 : & strten, targetpol)
564 : :
565 [ + - ]: 3 : ABI_SFREE(gred)
566 [ - + ]: 3 : ABI_FREE(d2flg)
567 : :
568 [ - + ]: 3 : end subroutine anaddb_driver_structural_response
569 : : !!***
570 : :
571 : : !!****f* m_anaddb_driver/anaddb_driver_susceptibilities
572 : : !! NAME
573 : : !! anaddb_driver_susceptibilities
574 : : !!
575 : : !! FUNCTION
576 : : !! Compute non-linear optical susceptibilities,
577 : : !! and if dtset%nlflag < 3, compute first-order change
578 : : !! in the linear dielectric susceptibility induced by an atomic displacement.
579 : : !! Then write susceptibilites to netcdf output.
580 : : !!
581 : : !! INPUTS
582 : : !!
583 : : !! OUTPUT
584 : : !!
585 : : !! SOURCE
586 : :
587 : 4 : subroutine anaddb_driver_susceptibilities(driver, dtset, ddb, ana_ncid, comm)
588 : :
589 : : !Arguments -------------------------------
590 : : class(anaddb_driver_type), intent(inout):: driver
591 : : type(anaddb_dataset_type), intent(in):: dtset
592 : : type(ddb_type), intent(in):: ddb
593 : : integer, intent(in):: ana_ncid
594 : : integer, intent(in):: comm
595 : :
596 : : !Local variables -------------------------------
597 : : integer, parameter:: master = 0
598 : : integer:: my_rank
599 : : integer:: ncerr
600 : : ! ************************************************************************
601 : :
602 [ - + ]: 4 : if (ddb%get_dchidet(dtset%ramansr, dtset%nlflag, driver%dchide, driver%dchidt) == 0) then
603 : 0 : ABI_ERROR("Cannot find block corresponding to non-linear optical susceptibilities in DDB file")
604 : : end if
605 : :
606 : : ! Save to the netcdf
607 : 4 : my_rank = xmpi_comm_rank(comm)
608 [ + - ]: 4 : if (my_rank == master) then
609 [ + + ]: 8 : ncerr = nctk_def_arrays(ana_ncid, [nctkarr_t("dchide", "dp", "three, three, three")], defmode=.True.)
610 [ - + ]: 4 : NCF_CHECK(ncerr)
611 [ - + ]: 4 : NCF_CHECK(nctk_set_datamode(ana_ncid))
612 [ - + ]: 4 : NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "dchide"), driver%dchide))
613 : :
614 : : ! dchidt only present if nlflag == 1 or 2
615 [ + + ]: 4 : if (dtset%nlflag < 3) then
616 : : ncerr = nctk_def_arrays(ana_ncid, [nctkarr_t("dchidt", "dp", &
617 [ + + ]: 6 : "number_of_atoms, three, three, three")], defmode=.True.)
618 [ - + ]: 3 : NCF_CHECK(ncerr)
619 [ - + ]: 3 : NCF_CHECK(nctk_set_datamode(ana_ncid))
620 [ - + ]: 3 : NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "dchidt"), driver%dchidt))
621 : : end if
622 : : end if
623 : :
624 : 4 : end subroutine anaddb_driver_susceptibilities
625 : : !!***
626 : :
627 : : !!****f* m_anaddb_driver/anaddb_driver_interatomic_force_constants
628 : : !! NAME
629 : : !! anaddb_driver_interatomic_force_constants
630 : : !!
631 : : !! FUNCTION
632 : : !! Interatomic forces calculation
633 : : !! Compute the interatomic force constants from a ddb.
634 : : !!
635 : : !! INPUTS
636 : : !!
637 : : !! OUTPUT
638 : : !!
639 : : !! SOURCE
640 : :
641 [ - + - + : 972 : subroutine anaddb_driver_interatomic_force_constants(driver, ifc, dtset, crystal, ddb, ana_ncid, comm)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + + + +
+ + + ]
642 : :
643 : : !Arguments -------------------------------
644 : : class(anaddb_driver_type), intent(inout):: driver
645 : : type(ifc_type), intent(out):: ifc
646 : : type(anaddb_dataset_type), intent(in):: dtset
647 : : type(crystal_t), intent(in):: crystal
648 : : type(ddb_type), intent(in):: ddb
649 : : integer, intent(in):: ana_ncid
650 : : integer, intent(in):: comm
651 : :
652 : : !Local variables -------------------------------
653 : : integer, parameter:: master = 0
654 : : integer:: my_rank
655 : : integer:: ii
656 [ + + + + : 864 : type(ifc_type):: Ifc_coarse
+ + - + ]
657 : : character(len = 500):: msg
658 : : integer:: ngqpt_coarse(3)
659 : : integer:: units(2)
660 : : integer:: qptrlatt(3,3), qptrlatt_coarse(3,3), nqshft, nqshft_coarse
661 : : real(dp):: q1shft(3,4), q1shft_coarse(3,4)
662 : : ! ************************************************************************
663 : :
664 [ + + ]: 162 : units = [std_out, ab_out]
665 : :
666 [ + + ]: 4428 : write(msg, '(a, a, (80a), a, a, a, a)' ) ch10, ('=',ii = 1, 80), ch10, ch10, &
667 [ + - ]: 4482 : ' Calculation of the interatomic forces ',ch10
668 : 54 : call wrtout(units, msg)
669 : :
670 : : ! Translate the (brav, ngqpt, q1shft) input into a brav-independent q-mesh lattice.
671 : 54 : call qmesh_from_brav_ngqpt(dtset%brav, dtset%ngqpt(1:3), dtset%nqshft, dtset%q1shft, qptrlatt, nqshft, q1shft)
672 : :
673 [ + + + + : 213 : if (any(dtset%qrefine(:) > 1)) then
+ + ]
674 : : ! Gaal-Nagy's algorithm in PRB 73 014117 [[cite:GaalNagy2006]]
675 : : ! Build the IFCs using the coarse q-mesh.
676 [ + + ]: 4 : do ii = 1, 3
677 : 4 : ngqpt_coarse(ii) = dtset%ngqpt(ii) / dtset%qrefine(ii)
678 : : end do
679 : : call qmesh_from_brav_ngqpt(dtset%brav, ngqpt_coarse, dtset%nqshft, dtset%q1shft, &
680 : 1 : qptrlatt_coarse, nqshft_coarse, q1shft_coarse)
681 : : call Ifc_coarse%init(crystal, ddb, &
682 : : dtset%asr, dtset%symdynmat, dtset%dipdip, dtset%rfmeth, qptrlatt_coarse, nqshft_coarse, q1shft_coarse, &
683 : : driver%epsinf, driver%zeff, driver%qdrp_cart, &
684 : : dtset%nsphere, dtset%rifcsph, dtset%prtsrlr, dtset%enunit, dtset%sys_dim, comm, &
685 : 1 : dipquad=dtset%dipquad, quadquad=dtset%quadquad, dielt_env=dtset%dielt_env,dielt_thick=dtset%dielt_thick)
686 : :
687 : : ! Now use the coarse q-mesh to fill the entries in dynmat(q)
688 : : ! on the dense q-mesh that cannot be obtained from the DDB file.
689 : : call ifc%init(crystal, ddb, &
690 : : dtset%asr, dtset%symdynmat, dtset%dipdip, dtset%rfmeth, &
691 : : qptrlatt, nqshft, q1shft, driver%epsinf, driver%zeff, driver%qdrp_cart, &
692 : : dtset%nsphere, dtset%rifcsph, dtset%prtsrlr, dtset%enunit, dtset%sys_dim, comm, &
693 : : Ifc_coarse=Ifc_coarse, dipquad=dtset%dipquad, quadquad=dtset%quadquad, &
694 : 1 : dielt_env=dtset%dielt_env, dielt_thick=dtset%dielt_thick)
695 : 1 : call Ifc_coarse%free()
696 : :
697 : : else
698 : : call ifc%init(crystal, ddb, &
699 : : dtset%asr, dtset%symdynmat, dtset%dipdip, dtset%rfmeth, &
700 : : qptrlatt, nqshft, q1shft, driver%epsinf, driver%zeff, driver%qdrp_cart, &
701 : : dtset%nsphere, dtset%rifcsph, dtset%prtsrlr, dtset%enunit, dtset%sys_dim, comm, &
702 : 53 : dipquad=dtset%dipquad, quadquad=dtset%quadquad, dielt_env=dtset%dielt_env, dielt_thick=dtset%dielt_thick)
703 : : end if
704 : :
705 [ + + ]: 108 : call ifc%print([std_out])
706 : :
707 : : ! If asr==6 (rotational invariance) we also need the estimation of the IFCs moments from real space.
708 [ + + - + ]: 54 : if (dtset%asr==6 .and. dtset%flexoflag==0) then
709 : 0 : write(msg, '(a, a)' )' Rotational invariance: will generate IFCs derivatives from IFCs moments',ch10
710 : 0 : call ifc%get_dcdq(crystal,driver%dcdq,driver%dcdqdq, comm)
711 : : ! If phi1 is not available from DDB, compute rotational invariance based on dcdq
712 : : ! For polar materials, we have to be extra cautious: indeed the typical electrostatics
713 : : ! model used to remove the non-analytical part is not invariant under rotation !
714 : : ! In 3D, the model leads to a divergent IFCs derivatives -> we cannot correct anything !
715 : : ! In 2D, the model leads to a finite torque and IFCs derivatives. We can add them back
716 : : ! during the imposition of rotational invariance...
717 : : end if
718 : :
719 : : ! Compute speed of sound.
720 [ - + ]: 54 : if (dtset%vs_qrad_tolkms(1) > zero) then
721 : 0 : call ifc%speedofsound(crystal, dtset%vs_qrad_tolkms, ana_ncid, comm)
722 : : end if
723 : :
724 : : ! Print analysis of the real-space interatomic force constants
725 : : ! TODO: ifc_out should not have side effects
726 : 54 : my_rank = xmpi_comm_rank(comm)
727 [ + - + + ]: 54 : if (my_rank == master .and. dtset%ifcout /= 0) then
728 : 13 : call ifc%write(dtset%ifcana, dtset%atifcflg, dtset%ifcout, dtset%prt_ifc, ana_ncid, dtset%prefix_outdata)
729 : : end if
730 : :
731 [ - + - + : 54 : end subroutine anaddb_driver_interatomic_force_constants
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
732 : : !!***
733 : :
734 : : !!****f* m_anaddb_driver/anaddb_driver_phdos
735 : : !! NAME
736 : : !! anaddb_driver_phdos
737 : : !!
738 : : !! FUNCTION
739 : : !! Compute phonon density of states.
740 : : !!
741 : : !! INPUTS
742 : : !!
743 : : !! OUTPUT
744 : : !!
745 : : !! SOURCE
746 : :
747 : 11 : subroutine anaddb_driver_phdos(driver, dtset, crystal, ifc, comm)
748 : :
749 : : !Arguments -------------------------------
750 : : class(anaddb_driver_type), intent(inout):: driver
751 : : type(anaddb_dataset_type), intent(in):: dtset
752 : : type(crystal_t), intent(in):: crystal
753 : : type(ifc_type), intent(in):: ifc
754 : : integer, intent(in):: comm
755 : :
756 : : !Local variables -------------------------------
757 : : integer, parameter:: master = 0
758 : : integer:: my_rank
759 : : integer:: ii
760 : : integer:: phdos_ncid, ncerr
761 : : character(len = fnlen):: phibz_prefix
762 : : character(len = 500):: msg
763 [ + + + + ]: 143 : type(phdos_t):: Phdos
764 : : integer:: units(2)
765 : : integer:: count_wminmax(2)
766 : : real(dp):: wminmax(2)
767 : : ! ************************************************************************
768 : :
769 : : ABI_UNUSED(driver%natom)
770 : 11 : my_rank = xmpi_comm_rank(comm)
771 [ + + ]: 33 : units = [std_out, ab_out]
772 : :
773 [ + - + + ]: 891 : write(msg, '(a, (80a), 4a)')ch10, ('=',ii = 1, 80), ch10, ch10, ' Calculation of phonon density of states ',ch10
774 : 11 : call wrtout(units, msg)
775 : :
776 : : ! Only 1 shift in q-mesh
777 : 11 : wminmax = zero
778 : 11 : phibz_prefix = trim(" ")
779 : : !phibz_prefix = "freq_displ" ! Uncomment this line to activate output of PHIBZ
780 : : ! ^^^^ GA: What the hell? FIXME ^^^^^
781 : 0 : do
782 : : call Phdos%init(crystal, Ifc, dtset%prtdos, dtset%dosdeltae, dtset%dossmear, dtset%ng2qpt, 1, dtset%q2shft, &
783 : 11 : phibz_prefix, wminmax, count_wminmax, comm, dos_maxmode=dtset%dos_maxmode)
784 [ + + + - : 33 : if (all(count_wminmax == 0)) exit
- + ]
785 : 0 : wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05; wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05
786 : 0 : call phdos%free()
787 : 0 : write(msg, "(a, 2f8.5)")"Initial frequency mesh not large enough. Recomputing PHDOS with wmin, wmax: ",wminmax
788 : 0 : call wrtout(std_out, msg)
789 : : end do
790 : :
791 [ + - ]: 11 : if (my_rank == master) then
792 : 11 : call phdos%print_msqd(dtset%prefix_outdata, dtset%ntemper, dtset%tempermin, dtset%temperinc)
793 : 11 : call phdos%print(strcat(dtset%prefix_outdata, "_PHDOS"))
794 : 11 : call phdos%print_debye(crystal%ucvol)
795 : 11 : call phdos%print_thermo(strcat(dtset%prefix_outdata, "_THERMO"), dtset%ntemper, dtset%tempermin, dtset%temperinc)
796 : :
797 : 11 : ncerr = nctk_open_create(phdos_ncid, strcat(dtset%prefix_outdata, "_PHDOS.nc"), xmpi_comm_self)
798 [ - + ]: 11 : NCF_CHECK_MSG(ncerr, "Creating PHDOS.nc file")
799 [ - + ]: 11 : NCF_CHECK(crystal%ncwrite(phdos_ncid))
800 : 11 : call phdos%ncwrite(phdos_ncid)
801 [ - + ]: 11 : NCF_CHECK(nf90_close(phdos_ncid))
802 : : end if
803 : :
804 : 11 : call phdos%free()
805 : :
806 [ - + - + : 11 : end subroutine anaddb_driver_phdos
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
807 : : !!***
808 : :
809 : : !!****f* m_anaddb_driver/anaddb_driver_thermal_supercell
810 : : !! NAME
811 : : !! anaddb_driver_thermal_supercell
812 : : !!
813 : : !! FUNCTION
814 : : !!
815 : : !! INPUTS
816 : : !!
817 : : !! OUTPUT
818 : : !!
819 : : !! SOURCE
820 : :
821 : 1 : subroutine anaddb_driver_thermal_supercell(driver, dtset, crystal, ifc)
822 : :
823 : : !Arguments -------------------------------
824 : : class(anaddb_driver_type), intent(in):: driver
825 : : type(anaddb_dataset_type), intent(in):: dtset
826 : : type(crystal_t), intent(in):: crystal
827 : : type(ifc_type), intent(in):: ifc
828 : :
829 : : !Local variables -------------------------------
830 : 1 : type(supercell_type), allocatable:: thm_scells(:)
831 : : ! ************************************************************************
832 : :
833 : : ABI_UNUSED(driver%natom)
834 [ - + + - : 13 : ABI_MALLOC(thm_scells, (dtset%ntemper))
+ - + - -
+ - + +
+ ]
835 : 1 : call zacharias_supercell_make(crystal, ifc, dtset%ntemper, dtset%thermal_supercell, dtset%tempermin, dtset%temperinc, thm_scells)
836 : 1 : call zacharias_supercell_print(dtset%prefix_outdata, dtset%ntemper, dtset%tempermin, dtset%temperinc, thm_scells)
837 : 1 : call thermal_supercell_free(dtset%ntemper, thm_scells)
838 [ + - + + : 11 : ABI_FREE(thm_scells)
- + - + -
+ - + - +
- + - + -
+ ]
839 : :
840 : 1 : end subroutine anaddb_driver_thermal_supercell
841 : : !!***
842 : :
843 : : !!****f* m_anaddb_driver/anaddb_driver_harmonic_thermo
844 : : !! NAME
845 : : !! anaddb_driver_harmonic_thermo
846 : : !!
847 : : !! FUNCTION
848 : : !! Phonon density of states and thermodynamical properties calculation
849 : : !!
850 : : !! INPUTS
851 : : !!
852 : : !! OUTPUT
853 : : !!
854 : : !! SOURCE
855 : :
856 : 4 : subroutine anaddb_driver_harmonic_thermo(driver, dtset, crystal, ifc, comm)
857 : :
858 : : !Arguments -------------------------------
859 : : class(anaddb_driver_type), intent(in):: driver
860 : : type(anaddb_dataset_type), intent(in):: dtset
861 : : type(crystal_t), intent(in):: crystal
862 : : type(ifc_type), intent(in):: ifc
863 : : integer, intent(in):: comm
864 : :
865 : : !Local variables -------------------------------
866 : : integer:: ii
867 : : character(len = 500):: msg
868 : : integer:: units(2)
869 : : ! ************************************************************************
870 : :
871 : : ABI_UNUSED(driver%natom)
872 [ + + ]: 12 : units = [std_out, ab_out]
873 : :
874 [ + + ]: 328 : write(msg, '(a, (80a), a, a, a, a, a, a, a, a)' ) ch10, ('=',ii = 1, 80), ch10, ch10, &
875 : 4 : ' Calculation of phonon density of states, ',ch10, &
876 : 4 : ' thermodynamical properties, ',ch10, &
877 [ + - ]: 332 : ' and Debye-Waller factors.',ch10
878 : 4 : call wrtout(units, msg)
879 : :
880 [ + - ]: 4 : if (dtset%thmflag == 1) then
881 : 4 : call harmonic_thermo(Ifc, crystal, ifc%amu, dtset, ab_out, dtset%prefix_outdata, comm)
882 : :
883 [ # # ]: 0 : else if (dtset%thmflag == 2) then
884 [ # # # # ]: 0 : write(msg, '(a, (80a), a, a, a, a)' ) ch10, ('=',ii = 1, 80), ch10, ch10, ' Entering thm9 routine with thmflag = 2 ',ch10
885 : 0 : call wrtout(units, msg)
886 : 0 : call harmonic_thermo(Ifc, crystal, ifc%amu, dtset, ab_out, dtset%prefix_outdata, comm, thmflag=dtset%thmflag)
887 : : end if
888 : :
889 : 4 : end subroutine anaddb_driver_harmonic_thermo
890 : : !!***
891 : :
892 : : !!****f* m_anaddb_driver/anaddb_driver_dielectric_q0
893 : : !! NAME
894 : : !! anaddb_driver_dielectric_q0
895 : : !!
896 : : !! FUNCTION
897 : : !! Compute dielectric tensor at Gamma and related properties:
898 : : !! mode effective charges, oscillator strength.
899 : : !!
900 : : !! INPUTS
901 : : !!
902 : : !! OUTPUT
903 : : !!
904 : : !! SOURCE
905 : :
906 : 27 : subroutine anaddb_driver_dielectric_q0(driver, dtset, crystal, ifc, ddb, asrq0, ana_ncid, comm)
907 : :
908 : : !Arguments -------------------------------
909 : : class(anaddb_driver_type), intent(inout):: driver
910 : : type(anaddb_dataset_type), intent(in):: dtset
911 : : type(crystal_t), intent(in):: crystal
912 : : type(ifc_type), intent(in):: ifc
913 : : type(ddb_type), intent(in):: ddb
914 : : type(asrq0_t), intent(inout):: asrq0
915 : : integer, intent(in):: ana_ncid
916 : : integer, intent(in):: comm
917 : :
918 : : !Local variables -------------------------------
919 : : integer:: ii, iblok
920 : : integer:: rfelfd(4), rfphon(4), rfstrs(4)
921 : : integer:: units(2)
922 : : real(dp) :: eta
923 : : character(len = 500):: msg
924 : : real(dp):: qphnrm(3), qphon(3, 3)
925 : 27 : real(dp), allocatable:: eigval(:,:)
926 : 27 : real(dp), allocatable:: eigvec(:,:,:,:,:)
927 : 27 : real(dp), allocatable:: lst(:)
928 : : ! ************************************************************************
929 : :
930 [ + + ]: 81 : units = [std_out, ab_out]
931 : :
932 [ - + + - : 81 : ABI_MALLOC(eigval, (3, driver%natom))
+ - + - -
+ - + ]
933 [ - + + - : 135 : ABI_MALLOC(eigvec, (2, 3, driver%natom, 3, driver%natom))
+ - - + +
- + - + -
- + - + ]
934 [ + - - + : 81 : ABI_MALLOC(lst, (dtset%nph2l+1))
+ - + - +
- - + -
+ ]
935 : :
936 [ + + ]: 91 : lst = zero
937 : :
938 : : !***************************************************************
939 : : ! Generates the dynamical matrix at Gamma
940 [ + + ]: 108 : qphon(:,1)=zero; qphnrm(1)=zero
941 : : ! Generation of the dynamical matrix in cartesian coordinates
942 [ + + ]: 27 : if (dtset%ifcflag == 1) then
943 : : ! Get d2cart using the interatomic forces and the
944 : : ! long-range coulomb interaction through Ewald summation
945 : : call gtdyn9(Ifc%atmfrc, driver%epsinf, Ifc%dipdip, &
946 : : Ifc%dyewq0, driver%d2cart, crystal%gmet, Ifc%gprimd, dtset%mpert, crystal%natom, &
947 : : Ifc%nrpt, qphnrm(1), qphon, crystal%rmet, Ifc%rprimd, Ifc%rpt, &
948 : : Ifc%trans, crystal%ucvol, Ifc%wghatm, crystal%xred, driver%zeff, driver%qdrp_cart, &
949 : : Ifc%ewald_option, eta, xmpi_comm_self, dtset%sys_dim,&
950 : : dipquad=Ifc%dipquad, quadquad=Ifc%quadquad, dielt_thick=dtset%dielt_thick,&
951 : 15 : dielt_env=dtset%dielt_env)
952 : :
953 [ + - ]: 12 : else if (dtset%ifcflag == 0) then
954 : : ! Look for the information in the DDB
955 [ + + + + : 108 : rfphon(1:2)=1; rfelfd(1:2)=2; rfstrs(1:2)=0
+ + ]
956 : 12 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, dtset%rfmeth)
957 [ + + ]: 12 : if (iblok == 0) then
958 [ + + + + ]: 3268 : driver%d2cart(:,1:dtset%msize)=zero
959 : : ! GA: I notice this situation happen in test tutorespfn[telast_3]
960 : : ! and I dont understand why the block is not found.
961 : : else
962 : : ! Copy the dynamical matrix in d2cart
963 [ + + + + ]: 898112 : driver%d2cart(:,1:dtset%msize)=ddb%val(:,:,iblok)
964 : : ! Eventually impose the acoustic sum rule
965 : 11 : call asrq0%apply(crystal%natom, dtset%mpert, dtset%msize, qphon, Crystal, driver%d2cart)
966 : : end if
967 : :
968 : : end if ! end of the generation of the dynamical matrix at gamma.
969 : : !***************************************************************
970 : :
971 : : ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
972 : : call dfpt_phfrq(ddb%amu, driver%displ, driver%d2cart, eigval, eigvec, crystal%indsym, &
973 : : dtset%mpert, crystal%nsym, crystal%natom, crystal%nsym, crystal%ntypat, driver%phfrq, qphnrm(1), qphon, &
974 : 27 : crystal%rprimd, dtset%symdynmat, crystal%symrel, crystal%symafm, crystal%typat, crystal%ucvol)
975 : :
976 : : ! calculation of the oscillator strengths, mode effective charge and
977 : : ! dielectric tensor, frequency dependent dielectric tensor (dieflag)
978 : : ! and mode by mode decomposition of epsilon if dieflag == 3
979 [ + + ]: 27 : if (dtset%dieflag /= 0) then
980 [ - + - - : 17 : if (driver%epsinf(1, 1)==one .and. driver%epsinf(2, 2)==one .and. driver%epsinf(3, 3)==one) then
- - ]
981 : 0 : write(msg, '(7a)') ch10, &
982 : 0 : ' The DDB file does not contain the derivatives w.r.t. electric field perturbation. ',ch10, &
983 : 0 : ' This is mandatory to calculate the dielectric constant, ',ch10, &
984 : 0 : ' Please check your DDB file or use dieflag = 0.',ch10
985 : 0 : ABI_ERROR(msg)
986 : : end if
987 : :
988 [ + - + + ]: 1377 : write(msg, '(a, (80a), a)' ) ch10, ('=',ii = 1, 80), ch10
989 : 17 : call wrtout(units, msg)
990 : :
991 : : ! Print the electronic contribution to the dielectric tensor
992 : : ! It can be extracted directly from the DDB if perturbation with E-field is present
993 : : call ddb_diel(crystal, ddb%amu, dtset, driver%dielt_rlx, driver%displ, driver%d2cart, driver%epsinf, driver%fact_oscstr, &
994 : 17 : ab_out, lst, dtset%mpert, crystal%natom, 0, driver%phfrq, comm, ana_ncid)
995 : : end if
996 : :
997 [ + - ]: 27 : ABI_SFREE(eigval)
998 [ + - ]: 27 : ABI_SFREE(eigvec)
999 [ + - ]: 27 : ABI_SFREE(lst)
1000 : :
1001 : 27 : end subroutine anaddb_driver_dielectric_q0
1002 : : !!***
1003 : :
1004 : : !!****f* m_anaddb_driver/anaddb_driver_nonlinear_response
1005 : : !! NAME
1006 : : !! anaddb_driver_nonlinear_response
1007 : : !!
1008 : : !! FUNCTION
1009 : : !! Non-linear response: electrooptic and Raman (q = Gamma, TO modes only)
1010 : : !!
1011 : : !! INPUTS
1012 : : !!
1013 : : !! OUTPUT
1014 : : !!
1015 : : !! SOURCE
1016 : :
1017 : 2 : subroutine anaddb_driver_nonlinear_response(driver, dtset, crystal, ana_ncid, comm)
1018 : :
1019 : : !Arguments -------------------------------
1020 : : class(anaddb_driver_type), intent(inout):: driver
1021 : : type(anaddb_dataset_type), intent(in):: dtset
1022 : : type(crystal_t), intent(in):: crystal
1023 : : integer, intent(in):: ana_ncid
1024 : : integer, intent(in):: comm
1025 : :
1026 : : !Local variables -------------------------------
1027 : : integer, parameter:: master = 0
1028 : : integer:: my_rank
1029 : : real(dp):: qphnrm(3), qphon(3, 3)
1030 : 2 : real(dp), allocatable:: rsus(:,:,:)
1031 : : ! ************************************************************************
1032 : :
1033 : 2 : my_rank = xmpi_comm_rank(comm)
1034 : :
1035 [ + - - + : 10 : ABI_MALLOC(rsus, (3*driver%natom, 3, 3))
+ - + - +
- + - + -
- + + - ]
1036 : :
1037 : : ! Raman susceptibilities for the 1st list (only TO modes at q = Gamma)
1038 [ + + ]: 8 : qphon(:,1)=zero
1039 : 2 : qphnrm(1)=zero
1040 : 2 : call ramansus(driver%d2cart, driver%dchide, driver%dchidt, driver%displ, dtset%mpert, crystal%natom, driver%phfrq, qphon, qphnrm(1), rsus, crystal%ucvol)
1041 : :
1042 [ + - ]: 2 : if (my_rank == master) then
1043 : 2 : call defwrite_raman_terms(ana_ncid, crystal%natom, rsus, driver%phfrq)
1044 : : end if
1045 : :
1046 : : ! EO coef:
1047 : 2 : call electrooptic(driver%dchide, dtset%dieflag, driver%epsinf, driver%fact_oscstr, crystal%natom, driver%phfrq, dtset%prtmbm, rsus, crystal%ucvol)
1048 : :
1049 [ + - ]: 2 : ABI_SFREE(rsus)
1050 : :
1051 : 2 : end subroutine anaddb_driver_nonlinear_response
1052 : : !!***
1053 : :
1054 : : !!****f* m_anaddb_driver/anaddb_driver_dielectric_nonana
1055 : : !! NAME
1056 : : !! anaddb_driver_dielectric_nonana
1057 : : !!
1058 : : !! FUNCTION
1059 : : !! Compute non-analyticity in the dielectric matrix and raman susceptibility
1060 : : !!
1061 : : !! INPUTS
1062 : : !!
1063 : : !! OUTPUT
1064 : : !!
1065 : : !! SOURCE
1066 : :
1067 : 24 : subroutine anaddb_driver_dielectric_nonana(driver, dtset, crystal, ddb, ana_ncid, comm)
1068 : :
1069 : : !Arguments -------------------------------
1070 : : class(anaddb_driver_type), intent(inout):: driver
1071 : : type(anaddb_dataset_type), intent(in):: dtset
1072 : : type(crystal_t), intent(in):: crystal
1073 : : type(ddb_type), intent(in):: ddb
1074 : : integer, intent(in):: ana_ncid
1075 : : integer, intent(in):: comm
1076 : :
1077 : : !Local variables -------------------------------
1078 : : integer, parameter:: master = 0
1079 : : integer:: ii, iphl2
1080 : : integer:: natom, nph2l
1081 : : integer:: my_rank
1082 : : integer:: units(2)
1083 : : character(len = 500):: msg
1084 : : real(dp):: qphnrm(3), qphon(3, 3)
1085 : 24 : real(dp), allocatable:: eigval(:,:)
1086 : 24 : real(dp), allocatable:: eigvec(:,:,:,:,:)
1087 : 24 : real(dp), allocatable:: rsus(:,:,:)
1088 : 24 : real(dp), allocatable:: lst(:)
1089 : : ! ************************************************************************
1090 : :
1091 : 24 : my_rank = xmpi_comm_rank(comm)
1092 [ + + ]: 72 : units = [std_out, ab_out]
1093 : :
1094 : 24 : natom = dtset%natom
1095 : 24 : nph2l = dtset%nph2l
1096 [ - + + - : 72 : ABI_MALLOC(eigval, (3, natom))
+ - + - -
+ - + ]
1097 [ - + + - : 120 : ABI_MALLOC(eigvec, (2, 3, natom, 3, natom))
+ - - + +
- + - + -
- + - + ]
1098 [ + - - + : 120 : ABI_MALLOC(rsus, (3*natom, 3, 3))
+ - + - +
- + - + -
- + - + ]
1099 [ + - - + : 72 : ABI_MALLOC(lst, (nph2l+1))
+ - + - +
- - + -
+ ]
1100 [ + + ]: 85 : lst = zero
1101 : :
1102 [ + - + + ]: 1944 : write(msg, '(a, (80a), a, a, a, a)' ) ch10, ('=',ii = 1, 80), ch10, ch10, ' Treat the second list of vectors ',ch10
1103 : 24 : call wrtout(units, msg)
1104 : :
1105 [ + - ]: 24 : if (my_rank == master) then
1106 : 24 : iphl2 = 0
1107 : 24 : call defwrite_nonana_terms(ana_ncid, iphl2, nph2l, dtset%qph2l, dtset%natom, driver%phfrq, driver%displ, "define")
1108 [ + + ]: 24 : if (dtset%nlflag == 1) then
1109 : 2 : call defwrite_nonana_raman_terms(ana_ncid, iphl2, nph2l, dtset%natom, rsus, "define")
1110 : : end if
1111 : : end if
1112 : :
1113 : : !Get the log of product of the square of the phonon frequencies without non-analyticities (q = 0)
1114 : : !For the Lyddane-Sachs-Teller relation, it is stored in lst(nph2+1)
1115 [ + + ]: 24 : if (dtset%dieflag /= 2 .and. dtset%dieflag /= 0) then
1116 [ + + ]: 108 : do ii = 4, 3*crystal%natom
1117 : 108 : lst(nph2l+1)=lst(nph2l+1)+2*log(driver%phfrq(ii))
1118 : : end do
1119 : : end if
1120 : :
1121 : : ! Examine every wavevector of this list
1122 [ + + ]: 61 : do iphl2 = 1, nph2l
1123 : :
1124 : : ! Initialisation of the phonon wavevector
1125 [ + + ]: 148 : qphon(:,1)=dtset%qph2l(:,iphl2)
1126 : 37 : qphnrm(1)=dtset%qnrml2(iphl2)
1127 : :
1128 : : ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
1129 : : ! for the second list of wv (can include non-analyticities if q /= 0)
1130 : : call dfpt_phfrq(ddb%amu, driver%displ, driver%d2cart, eigval, eigvec, crystal%indsym, &
1131 : : dtset%mpert, crystal%nsym, crystal%natom, crystal%nsym, crystal%ntypat, driver%phfrq, qphnrm(1), qphon, &
1132 : 37 : crystal%rprimd, dtset%symdynmat, crystal%symrel, crystal%symafm, crystal%typat, crystal%ucvol)
1133 : :
1134 : : ! Write the phonon frequencies for the second list of wv (can include non-analyticities if q /= 0)
1135 : 37 : call dfpt_prtph(driver%displ, dtset%eivec, dtset%enunit, ab_out, dtset%natom, driver%phfrq, qphnrm(1), qphon)
1136 : : ! TODO: Mode effective charge could be printed here for LO modes (EB)
1137 : :
1138 [ + - ]: 37 : if (my_rank == master) then
1139 : : ! Loop is not MPI-parallelized--> no need for MPI-IO API.
1140 : 37 : call defwrite_nonana_terms(ana_ncid, iphl2, nph2l, dtset%qph2l, dtset%natom, driver%phfrq, driver%displ, "write")
1141 : : end if
1142 : :
1143 : : ! Get the log of product of the square of the phonon frequencies with non-analyticities (q-->0)
1144 : : ! for the Lyddane-Sachs-Teller relation
1145 : : ! The fourth mode should have positive frequency otherwise there is an instability: LST relationship should not be evaluated
1146 : : ! Isn't it tested somewhere else (i.e. stop of the code if there are imaginary freq.)? (EB)
1147 [ + + ]: 37 : if (dtset%dieflag /= 2 .and. dtset%dieflag /= 0) then
1148 [ + + ]: 149 : do ii = 4, 3*crystal%natom
1149 : 149 : lst(iphl2)=lst(iphl2)+2*log(driver%phfrq(ii))
1150 : : end do
1151 : : end if
1152 : :
1153 : : ! Write Raman susceptibilities for the 2nd list (can includes LO modes if q /= 0 0 0)
1154 [ + + ]: 61 : if (dtset%nlflag == 1) then
1155 : 4 : call ramansus(driver%d2cart, driver%dchide, driver%dchidt, driver%displ, dtset%mpert, crystal%natom, driver%phfrq, qphon, qphnrm(1), rsus, crystal%ucvol)
1156 [ + - ]: 4 : if (my_rank == master) then
1157 : 4 : call defwrite_nonana_raman_terms(ana_ncid, iphl2, nph2l, dtset%natom, rsus, "write")
1158 : : end if
1159 : : end if ! nlflag = 1 (Raman suscep for the 2nd list of wv.)
1160 : : end do ! iphl2
1161 : :
1162 : : ! Lyddane-Sachs-Teller relation:
1163 [ + + ]: 24 : if (dtset%dieflag /= 2 .and. dtset%dieflag /= 0) then
1164 : : call ddb_diel(crystal, ddb%amu, dtset, driver%dielt_rlx, driver%displ, driver%d2cart, driver%epsinf, driver%fact_oscstr, &
1165 : 9 : ab_out, lst, dtset%mpert, crystal%natom, nph2l, driver%phfrq, comm, ana_ncid)
1166 : : end if
1167 : :
1168 [ + - ]: 24 : ABI_SFREE(eigval)
1169 [ + - ]: 24 : ABI_SFREE(eigvec)
1170 [ + - ]: 24 : ABI_SFREE(rsus)
1171 [ + - ]: 24 : ABI_SFREE(lst)
1172 : :
1173 : 24 : end subroutine anaddb_driver_dielectric_nonana
1174 : : !!***
1175 : :
1176 : : !!****f* m_anaddb_driver/anaddb_driver_internal_strain
1177 : : !! NAME
1178 : : !! anaddb_driver_internal_strain
1179 : : !!
1180 : : !! FUNCTION
1181 : : !!
1182 : : !! INPUTS
1183 : : !!
1184 : : !! OUTPUT
1185 : : !!
1186 : : !! SOURCE
1187 : :
1188 : 6 : subroutine anaddb_driver_internal_strain(driver, dtset, ddb, asrq0)
1189 : :
1190 : : !Arguments -------------------------------
1191 : : class(anaddb_driver_type), intent(inout):: driver
1192 : : type(anaddb_dataset_type), intent(in):: dtset
1193 : : type(ddb_type), intent(in):: ddb
1194 : : type(asrq0_t), intent(in):: asrq0
1195 : :
1196 : : !Local variables -------------------------------
1197 : : integer:: ii, iblok
1198 : : integer:: prt_internalstr
1199 : : integer:: units(2)
1200 : : integer:: rfelfd(4), rfphon(4), rfstrs(4)
1201 : : character(len = 500):: msg
1202 : : real(dp):: qphnrm(3), qphon(3, 3)
1203 : : ! ************************************************************************
1204 : :
1205 [ + + ]: 18 : units = [std_out, ab_out]
1206 : :
1207 : : ! Here treating the internal strain tensors at Gamma point
1208 [ + + ]: 492 : write(msg, '(a, a, (80a), a, a, a, a)') ch10, ('=',ii = 1, 80), ch10, ch10, &
1209 [ + - ]: 498 : ' Calculation of the internal-strain tensor',ch10
1210 : 6 : call wrtout(units, msg)
1211 : :
1212 [ + - ]: 6 : if (dtset%instrflag == 1) then
1213 : 6 : call wrtout(std_out, 'instrflag = 1, so extract the internal strain constant from the 2DTE')
1214 : :
1215 : : ! looking after the no. of blok that contains the internal strain tensor
1216 [ + + ]: 24 : qphon(:,1)=zero; qphnrm(1)=zero
1217 [ + + + + : 54 : rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
+ + ]
1218 : :
1219 : 6 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, dtset%rfmeth)
1220 [ - + ]: 6 : if (iblok == 0) then
1221 : 0 : ABI_ERROR("DDB file must contain both uniaxial and shear strain for piezoelectric, Check your calculations")
1222 : : end if
1223 : :
1224 : : ! then print the internal stain tensor
1225 : 6 : prt_internalstr = 2
1226 : 6 : call ddb_internalstr(dtset%asr, ddb%val, asrq0%d2asr, iblok, driver%instrain, ab_out, dtset%mpert, ddb%natom, ddb%nblok, prt_internalstr)
1227 : : end if
1228 : :
1229 : 6 : end subroutine anaddb_driver_internal_strain
1230 : : !!***
1231 : :
1232 : : !!****f* m_anaddb_driver/anaddb_driver_elastic_tensor
1233 : : !! NAME
1234 : : !! anaddb_driver_elastic_tensor
1235 : : !!
1236 : : !! FUNCTION
1237 : : !!
1238 : : !! INPUTS
1239 : : !!
1240 : : !! OUTPUT
1241 : : !!
1242 : : !! SOURCE
1243 : :
1244 : 6 : subroutine anaddb_driver_elastic_tensor(driver, dtset, crystal, ddb, asrq0, ana_ncid)
1245 : :
1246 : : !Arguments -------------------------------
1247 : : class(anaddb_driver_type), intent(inout):: driver
1248 : : type(anaddb_dataset_type), intent(in):: dtset
1249 : : type(crystal_t), intent(in):: crystal
1250 : : type(ddb_type), intent(in):: ddb
1251 : : type(asrq0_t), intent(inout):: asrq0
1252 : : integer, intent(in):: ana_ncid
1253 : :
1254 : : !Local variables -------------------------------
1255 : : integer:: ii, iblok, iblok_stress
1256 : : integer:: units(2)
1257 : : character(len = 500):: msg
1258 : : integer:: rfelfd(4), rfphon(4), rfstrs(4)
1259 : : real(dp):: qphnrm(3), qphon(3, 3)
1260 : : real(dp):: compl(6, 6), compl_clamped(6, 6), compl_stress(6, 6)
1261 : : real(dp):: elast_clamped(6, 6), elast_stress(6, 6)
1262 : : ! ************************************************************************
1263 : :
1264 [ + + ]: 18 : units = [std_out, ab_out]
1265 : :
1266 : : ! here treating the elastic tensors at Gamma Point
1267 [ + + ]: 492 : write(msg, '(a, a, (80a), a, a, a, a, a, a)') ch10, ('=',ii = 1, 80), ch10, ch10, &
1268 [ + - ]: 498 : ' Calculation of the elastic and compliances tensor (Voigt notation)',ch10
1269 : 6 : call wrtout(units, msg)
1270 : :
1271 : 6 : call wrtout(std_out, 'so extract the elastic constant from the 2DTE')
1272 : :
1273 : : ! look after the blok no. that contains the stress tensor
1274 [ + + ]: 24 : qphon(:,1)=zero; qphnrm(1)=zero
1275 [ + + + + : 54 : rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=0
+ + ]
1276 : :
1277 : 6 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, BLKTYP_d1E_xx)
1278 : 6 : iblok_stress = iblok
1279 : :
1280 : : ! look after the blok no.iblok that contains the elastic tensor
1281 [ + + ]: 24 : qphon(:,1)=zero; qphnrm(1)=zero
1282 [ + + + + : 54 : rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
+ + ]
1283 : :
1284 : : ! for both diagonal and shear parts
1285 : 6 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, dtset%rfmeth)
1286 [ - + ]: 6 : if (iblok == 0) then
1287 : 0 : ABI_ERROR("DDB file must contain both uniaxial and shear strain when elaflag != 0, Check your calculations")
1288 : : end if
1289 : :
1290 : : ! print the elastic tensor
1291 : : call ddb_elast(dtset, crystal, ddb%val, compl, compl_clamped, compl_stress, asrq0%d2asr, &
1292 : : driver%elast, elast_clamped, elast_stress, iblok, iblok_stress, &
1293 : 6 : driver%instrain, ab_out, dtset%mpert, crystal%natom, ddb%nblok, ana_ncid)
1294 : :
1295 : 6 : end subroutine anaddb_driver_elastic_tensor
1296 : : !!***
1297 : :
1298 : : !!****f* m_anaddb_driver/anaddb_driver_piezoelectric_tensor
1299 : : !! NAME
1300 : : !! anaddb_driver_piezoelectric_tensor
1301 : : !!
1302 : : !! FUNCTION
1303 : : !!
1304 : : !! INPUTS
1305 : : !!
1306 : : !! OUTPUT
1307 : : !!
1308 : : !! SOURCE
1309 : :
1310 : 5 : subroutine anaddb_driver_piezoelectric_tensor(driver, dtset, crystal, ddb, ana_ncid)
1311 : :
1312 : : !Arguments -------------------------------
1313 : : class(anaddb_driver_type), intent(inout):: driver
1314 : : type(anaddb_dataset_type), intent(in):: dtset
1315 : : type(crystal_t), intent(in):: crystal
1316 : : type(ddb_type), intent(in):: ddb
1317 : : integer, intent(in):: ana_ncid
1318 : :
1319 : : !Local variables -------------------------------
1320 : : integer:: ii, iblok
1321 : : integer:: units(2)
1322 : : character(len = 500):: msg
1323 : : integer:: rfelfd(4), rfphon(4), rfstrs(4)
1324 : : real(dp):: qphnrm(3), qphon(3, 3)
1325 : : real(dp):: piezo(6, 3)
1326 : : ! ************************************************************************
1327 : :
1328 [ + + ]: 15 : units = [std_out, ab_out]
1329 : :
1330 : : ! Here treating the piezoelectric tensor at Gamma Point
1331 [ + + ]: 410 : write(msg, '(a, a, (80a), a, a, a, a, a)') ch10, ('=',ii = 1, 80), ch10, ch10, &
1332 : 5 : ' Calculation of the tensor related to piezoelectric effetc',ch10, &
1333 [ + - ]: 415 : ' (Elastic indices in Voigt notation)',ch10
1334 : 5 : call wrtout(units, msg)
1335 : :
1336 : 5 : call wrtout(std_out, 'extract the piezoelectric constant from the 2DTE')
1337 : :
1338 : : ! Looking for the gamma point block
1339 [ + + ]: 20 : qphon(:,1)=zero; qphnrm(1)=zero
1340 [ + + + + : 45 : rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
+ + ]
1341 : :
1342 : : ! For both diagonal and shear parts
1343 : 5 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, dtset%rfmeth)
1344 [ - + ]: 5 : if (iblok == 0) then
1345 : 0 : ABI_ERROR("DDB file must contain both uniaxial and shear strain for piezoelectric, Check your calculations")
1346 : : end if
1347 : :
1348 : : ! Then print out the piezoelectric constants
1349 : : call ddb_piezo(dtset, ddb%val, driver%dielt_rlx, driver%elast, iblok, &
1350 : : & driver%instrain, ab_out, dtset%mpert, crystal%natom, ddb%nblok, piezo, &
1351 : 5 : & crystal%ucvol, ana_ncid)
1352 : :
1353 : 5 : end subroutine anaddb_driver_piezoelectric_tensor
1354 : : !!***
1355 : :
1356 : : !!****f* m_anaddb_driver/anaddb_driver_flexoelectric_tensor
1357 : : !! NAME
1358 : : !! anaddb_driver_flexoelectric_tensor
1359 : : !!
1360 : : !! FUNCTION
1361 : : !!
1362 : : !! INPUTS
1363 : : !!
1364 : : !! OUTPUT
1365 : : !!
1366 : : !! SOURCE
1367 : :
1368 : 2 : subroutine anaddb_driver_flexoelectric_tensor(driver, dtset, crystal, ddb, ddb_lw, ddb_hdr, asrq0)
1369 : :
1370 : : !Arguments -------------------------------
1371 : : class(anaddb_driver_type), intent(in):: driver
1372 : : type(anaddb_dataset_type), intent(in):: dtset
1373 : : type(crystal_t), intent(in):: crystal
1374 : : type(ddb_type), intent(in):: ddb, ddb_lw
1375 : : type(ddb_hdr_type), intent(in):: ddb_hdr
1376 : : type(asrq0_t), intent(in):: asrq0
1377 : :
1378 : : !Local variables -------------------------------
1379 : : integer:: ii
1380 : : integer:: units(2)
1381 : : character(len = 500):: msg
1382 : : ! ************************************************************************
1383 : :
1384 [ + + ]: 6 : units = [std_out, ab_out]
1385 : :
1386 : : ! Here treating the flexoelectric tensor
1387 [ + + ]: 164 : write(msg, '(a, a, (80a), a, a, a, a)') ch10, ('=',ii = 1, 80), ch10, ch10, &
1388 [ + - ]: 166 : ' Calculation of the tensors related to flexoelectric effect',ch10
1389 : 2 : call wrtout(units, msg)
1390 : :
1391 : : ! Compute and print the contributions to the flexoelectric tensor
1392 : : call ddb_flexo(dtset%asr, asrq0%d2asr, ddb, ddb_lw, ddb_hdr%ddb_version, crystal, &
1393 : 2 : & dtset%filename_ddb, dtset%flexoflag, dtset%prtvol, driver%zeff)
1394 : :
1395 : 2 : end subroutine anaddb_driver_flexoelectric_tensor
1396 : :
1397 : : !!***
1398 : :
1399 : : !!****f* m_anaddb_driver/anaddb_driver_get_dcdq
1400 : : !! NAME
1401 : : !! anaddb_driver_get_dcdq
1402 : : !!
1403 : : !! FUNCTION
1404 : : !!
1405 : : !! INPUTS
1406 : : !!
1407 : : !! OUTPUT
1408 : : !!
1409 : : !! SOURCE
1410 : :
1411 : 1 : subroutine anaddb_driver_get_dcdq(driver, dtset, ddb, ddb_lw, ddb_hdr)
1412 : :
1413 : : !Arguments -------------------------------
1414 : : class(anaddb_driver_type), intent(inout):: driver
1415 : : type(anaddb_dataset_type), intent(in):: dtset
1416 : : type(ddb_type), intent(in):: ddb, ddb_lw
1417 : : type(ddb_hdr_type), intent(in):: ddb_hdr
1418 : :
1419 : : !Local variables -------------------------------
1420 : : integer:: ii
1421 : : integer:: units(2)
1422 : : character(len = 500):: msg
1423 : : !real(dp):: dcdq(3,ddb%natom,3,ddb%natom,3)
1424 : :
1425 : : ! ************************************************************************
1426 : :
1427 [ + + ]: 3 : units = [std_out, ab_out]
1428 : :
1429 : : ! Here treating the flexoelectric tensor
1430 [ + + ]: 82 : write(msg, '(a, a, (80a), a, a, a, a)') ch10, ('=',ii = 1, 80), ch10, ch10, &
1431 [ + - ]: 83 : ' Extraction from DDB of the IFCs derivatives',ch10
1432 : 1 : call wrtout(units, msg)
1433 : :
1434 : : ! Compute and print the contributions to the flexoelectric tensor
1435 : 1 : call ddb_phi1(ddb,ddb_lw,ddb_hdr%ddb_version, dtset%filename_ddb, driver%dcdq, driver%dcdqdq,ddb%natom)
1436 : :
1437 : 1 : end subroutine anaddb_driver_get_dcdq
1438 : :
1439 : : !!***
1440 : :
1441 : : !!****f* m_anaddb_driver/anaddb_driver_lattice_wannier
1442 : : !! NAME
1443 : : !! anaddb_driver_lattice_wannier
1444 : : !!
1445 : : !! FUNCTION
1446 : : !! Lattice Wannier function calculation.
1447 : : !! Compute the Dynamical matrix for a dense Q-mesh
1448 : : !! Input the eigenvectors and eigenvalues to the Lattcie Wannier module
1449 : : !! Construct the Lattice Wannier functions
1450 : : !!
1451 : : !! INPUTS
1452 : : !!
1453 : : !! OUTPUT
1454 : : !!
1455 : : !! SOURCE
1456 : :
1457 : 3 : subroutine anaddb_driver_lattice_wannier(driver, dtset, crystal, ifc, comm)
1458 : :
1459 : : !Arguments -------------------------------
1460 : : class(anaddb_driver_type), intent(in):: driver
1461 : : type(anaddb_dataset_type), intent(in):: dtset
1462 : : type(crystal_t),intent(in):: crystal
1463 : : type(ifc_type), intent(in):: ifc
1464 : : integer,intent(in):: comm
1465 : :
1466 : : !Local variables -------------------------------
1467 : : integer:: ii
1468 : : character(len = 500):: msg
1469 : : integer:: units(2)
1470 : : ! ************************************************************************
1471 : :
1472 [ + + ]: 9 : units = [std_out, ab_out]
1473 : : ABI_UNUSED(driver%natom)
1474 : :
1475 [ + - + + ]: 243 : write(msg, '(a, (80a), 4a)')ch10, ('=',ii = 1, 80), ch10, ch10, ' Calculation of lattice Wannier functions ',ch10
1476 : 3 : call wrtout(units, msg)
1477 : 3 : call run_lattice_wannier(ifc=ifc, crystal=crystal, dtset=dtset, prefix=dtset%prefix_outdata, comm=comm)
1478 [ + - + + ]: 243 : write(msg, '(a, (80a))')ch10, ('=',ii = 1, 80)
1479 : 3 : call wrtout(units, msg)
1480 : :
1481 : 3 : end subroutine anaddb_driver_lattice_wannier
1482 : : !!***
1483 : :
1484 : : !!***
1485 : :
1486 : : !!****f* m_anaddb_driver/anaddb_driver_convertdim_dielt
1487 : : !! NAME
1488 : : !! anaddb_driver_convertdim_dielt
1489 : : !!
1490 : : !! FUNCTION
1491 : : !! Supposing a reduced dimensionality of the problem, convert DDB dielectric tensor
1492 : : !! to the low-dimensionality dielectric tensor of the isolated material. Indeed,
1493 : : !! first-principles calculations requires vacuum buffer(s) to isolate the low-dimensional
1494 : : !! materials, that spuriously contributes to the dielectric response of the whole slab
1495 : : !!
1496 : : !! INPUTS
1497 : : !!
1498 : : !! OUTPUT
1499 : : !!
1500 : : !! SOURCE
1501 : :
1502 : 2 : subroutine anaddb_driver_convertdim_dielt(driver, rprimd, sys_dim, dielt_thick)
1503 : :
1504 : : !Arguments -------------------------------
1505 : : class(anaddb_driver_type), intent(inout):: driver
1506 : : integer,intent(in):: sys_dim
1507 : : real(dp) :: dielt_thick(2), rprimd(3,3)
1508 : :
1509 : : !Local variables -------------------------------
1510 : : integer:: idir, jdir, kdir, bool_isolated(3)
1511 : : real(dp) :: thick, epsinf_conv(3,3), acell(3)
1512 : : character(len = 500):: msg
1513 : :
1514 : : ! ************************************************************************
1515 : 2 : epsinf_conv=zero
1516 : 2 : thick = dielt_thick(1)-dielt_thick(2)
1517 : 2 : bool_isolated(:) = 0
1518 : : ! When periodic, additional variable spaces coming from dynamical matrices derivatives
1519 [ + - ]: 2 : if (sys_dim == 1) then ! 3D
1520 : : bool_isolated = 0
1521 [ - + ]: 2 : elseif (sys_dim == 2) then ! 2D yz
1522 : 0 : bool_isolated(1) = 1
1523 [ - + ]: 2 : elseif (sys_dim == 3) then ! 2D xz
1524 : 0 : bool_isolated(2) = 1
1525 [ + - ]: 2 : elseif (sys_dim == 4) then ! 2D xy
1526 : 2 : bool_isolated(3) = 1
1527 : : else
1528 : : write(msg,'(3a,i0)') &
1529 : 0 : 'For dipole-dipole in 2D, the argument sys_dim should',ch10,&
1530 : 0 : 'be between 1 and 4. However, sys_dim = ',sys_dim
1531 : 0 : ABI_ERROR(msg)
1532 : : end if
1533 [ + + ]: 8 : do idir = 1,3
1534 [ + + ]: 26 : do jdir =1,3
1535 [ + + + + : 24 : if (bool_isolated(idir)==1 .and. bool_isolated(jdir)==0 .and. driver%epsinf(idir,jdir)>tol2) then
- + ]
1536 : : write(msg,'(3a)')&
1537 : 0 : 'Along a confined + periodic directions, the dielectric tensor should be 0.',ch10,&
1538 : 0 : 'However, a component is found to be larger than 1e-6. Please check your DDB file'
1539 : 0 : ABI_ERROR(msg)
1540 : : end if
1541 : : end do
1542 : : end do
1543 [ + + ]: 8 : do idir=1,3
1544 [ + + ]: 26 : acell(idir)= sqrt(dot_product(rprimd(idir,:),rprimd(idir,:)))
1545 : : end do
1546 : : write(msg,'(2a)') &
1547 : 2 : 'Conversion of the dielectric tensor for 2D materials to account for vacuum',ch10
1548 [ + + ]: 6 : call wrtout([std_out, ab_out],msg)
1549 [ + + ]: 8 : do kdir=1,3
1550 [ + + ]: 8 : if (bool_isolated(kdir)==1) then
1551 [ + + ]: 8 : do idir=1,3
1552 [ + + ]: 26 : do jdir=1,3
1553 [ + + + + ]: 24 : if (bool_isolated(idir)==0 .and. bool_isolated(jdir)==0) then
1554 : : ! Capacitors in parallel for periodic directions between 2D and vacuum
1555 [ + + ]: 8 : if (idir==jdir) then
1556 : 4 : epsinf_conv(idir,jdir)=one-acell(kdir)/dielt_thick(1)*(one-driver%epsinf(idir,jdir))
1557 : : else
1558 : 4 : epsinf_conv(idir,jdir)=acell(kdir)/dielt_thick(1)*driver%epsinf(idir,jdir)
1559 : : end if
1560 [ + + + + ]: 10 : elseif (bool_isolated(idir)==1 .and. bool_isolated(jdir)==1) then
1561 : : ! Capacitors in series for confined direction between 2D and vacuum
1562 : 2 : epsinf_conv(idir,jdir)= one/(one-acell(kdir)/thick*(one-one/driver%epsinf(idir,jdir)))
1563 : : end if
1564 : : end do
1565 : : end do
1566 : : end if
1567 : : end do
1568 : : write(msg,'(a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a)')&
1569 : 2 : 'Converted dielectric tensor',epsinf_conv(1,1),' ',epsinf_conv(1,2), ' ', epsinf_conv(1,3),' '&
1570 : 2 : ,epsinf_conv(2,1),' ',epsinf_conv(2,2), ' ', epsinf_conv(2,3),' '&
1571 : 4 : ,epsinf_conv(3,1),' ',epsinf_conv(3,2), ' ', epsinf_conv(3,3), ch10
1572 [ + + ]: 6 : call wrtout([std_out, ab_out],msg)
1573 [ + + + + ]: 26 : driver%epsinf = epsinf_conv
1574 : :
1575 : 2 : end subroutine anaddb_driver_convertdim_dielt
1576 : :
1577 [ # # # # : 0 : end module m_anaddb_driver
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1578 : : !!***
|
