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