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 : !!***
|
