Line data Source code
1 : !!****m* ABINIT/m_ddb_diel
2 : !! NAME
3 : !! m_ddb_diel
4 : !!
5 : !! FUNCTION
6 : !! This module provides routines for the calculation of the dielectric constant (anaddb)
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XG,XW, MVeithen, EB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_ddb_diel
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_xmpi
27 : use m_abicore
28 : use m_ddb
29 : use m_nctk
30 : use netcdf
31 :
32 : use m_anaddb_dataset, only : anaddb_dataset_type
33 : use m_crystal, only : crystal_t
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : public :: ddb_diel
41 : public :: alignph
42 : !!***
43 :
44 : contains
45 : !!***
46 :
47 : !!****f* ABINIT/ddb_diel
48 : !!
49 : !! NAME
50 : !! ddb_diel
51 : !!
52 : !! FUNCTION
53 : !! Get the frequency-dependent dielectric matrix, as well as the
54 : !! oscillator strengths and mode effective charges,
55 : !! and reflectivities (without damping)
56 : !! See the definitions Eq.(53-54) in PRB55, 10355 (1997) [[cite:Gonze1997a]].
57 : !!
58 : !! INPUTS
59 : !! amu(ntypat)=mass of the atoms (atomic mass unit)
60 : !! anaddb_dtset= (derived datatype) contains all the input variables
61 : !! matrix (diagonal in the atoms)
62 : !! displ(2,3*natom,3*natom)=
63 : !! the displacements of atoms in cartesian coordinates.
64 : !! The first index means either the real or the imaginary part,
65 : !! The second index runs on the direction and the atoms displaced
66 : !! The third index runs on the modes.
67 : !! d2cart(2,3,mpert,3,mpert)=dynamical matrix, effective charges, dielectric tensor,... all in cartesian coordinates
68 : !! iout=unit number for outputs
69 : !! lst(3*nph2l)=log. of product of frequencies**2, needed to calculate
70 : !! the generalized Lyddane-Sachs-Teller relation at zero frequency
71 : !! mpert =maximum number of ipert
72 : !! natom=number of atoms in unit cell
73 : !! nph2l=input variable from anaddb_dtset, needed to dimension lst
74 : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
75 : !! matrix eigenvalues, except if these are negative, and in this
76 : !! case, give minus the square root of the absolute value
77 : !! of the matrix eigenvalues). Hartree units.
78 : !! comm=MPI communicator.
79 : !! ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
80 : !!
81 : !! OUTPUT
82 : !! fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes,
83 : !! for different direction of the electric field;
84 : !! dielt_rlx(3,3) relaxed ion (zero frequency) dielectric tensor.
85 : !!
86 : !! NOTES
87 : !! 1. The phonon frequencies phfrq should correspond to the
88 : !! wavevector at Gamma, without any non-analyticities.
89 : !! 2. Should clean for no imaginary part ...
90 : !! This routine should be used only by one processor.
91 : !! 3. frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
92 : !! mode effective charges for the different eigenmodes,
93 : !! for different direction of the electric field
94 : !!
95 : !! SOURCE
96 :
97 26 : subroutine ddb_diel(Crystal,amu,anaddb_dtset,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
98 26 : & iout,lst,mpert,natom,nph2l,phfrq,comm,ncid)
99 :
100 : !Arguments -------------------------------
101 : !scalars
102 : integer,intent(in) :: iout,mpert,natom,nph2l,comm,ncid
103 : type(crystal_t),intent(in) :: Crystal
104 : type(anaddb_dataset_type),intent(in) :: anaddb_dtset
105 : real(dp),intent(in) :: lst(nph2l+1)
106 :
107 : !arrays
108 : real(dp),intent(in) :: amu(Crystal%ntypat),d2cart(2,3,mpert,3,mpert)
109 : real(dp),intent(in) :: phfrq(3*natom),epsinf(3,3)
110 : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
111 : real(dp),intent(out) :: dielt_rlx(3,3),fact_oscstr(2,3,3*natom)
112 :
113 : !Local variables -------------------------
114 : !scalars
115 : integer,parameter :: master=0
116 : integer :: dieflag,idir1,idir2,ifreq,ii,imode,iphl2,nfreq
117 : integer :: nprocs,my_rank,ncerr
118 : real(dp) :: afreq,difffr,eps,q2,ucvol
119 : character(len=500) :: message
120 : !arrays
121 : real(dp) :: qphon(3),refl(3)
122 26 : real(dp),allocatable :: frdiel(:,:,:),modez(:,:,:),oscstr(:,:,:,:),dielt_modedecompo(:,:,:)
123 :
124 : ! *********************************************************************
125 :
126 26 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
127 :
128 26 : dieflag=anaddb_dtset%dieflag
129 26 : nfreq=anaddb_dtset%nfreq
130 :
131 26 : ucvol = Crystal%ucvol
132 :
133 :
134 : ! frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
135 : ! modez=mode effective charge
136 78 : ABI_MALLOC(frdiel,(3,3,nfreq))
137 104 : ABI_MALLOC(oscstr,(2,3,3,3*natom))
138 78 : ABI_MALLOC(modez,(2,3,3*natom))
139 78 : ABI_MALLOC(dielt_modedecompo,(3,3,3*natom))
140 :
141 : !In case only the electronic dielectric constant should be printed
142 26 : if (dieflag==2) then
143 7 : call ddb_diel_elec(iout,epsinf)
144 : else
145 : ! In case the ionic contribution to the dielectric tensor is asked
146 19 : if (dieflag/=2 .and. nph2l==0) then
147 : ! Check if the alignement of phonon modes eigenvector is requested from the input flag alphon;
148 : ! useful in case of degenerate modes
149 10 : if (anaddb_dtset%alphon > 0) then
150 : write(message, '(3a)' )&
151 3 : ' The alphon input variable is non-zero, will mix the degenerate phonon modes',ch10,&
152 6 : ' in order to align the mode effective charges with the cartesian axes.'
153 3 : call wrtout(std_out,message,'COLL')
154 3 : call wrtout(iout,message,'COLL')
155 3 : call alignph(amu,displ,d2cart,mpert,natom,Crystal%ntypat,phfrq,Crystal%typat)
156 : end if ! alignment of the phonon eigenvectors
157 :
158 : ! Compute the mode effective charge and oscillator strength
159 10 : call ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank)
160 :
161 : ! Calculation of epsilon_r (Eq.55 PRB 55, 10355)
162 : ! Check the acousticity of the three lowest modes, assuming
163 : ! that they are ordered correctly
164 10 : if (abs(phfrq(1))>abs(phfrq(4)))then
165 : ! This means that there is at least one mode with truly negative frequency
166 : write(message, '(12a,4es16.8)' )&
167 0 : 'The lowest mode appears to be a "true" negative mode,',ch10,&
168 0 : 'and not an acoustic mode. This precludes the computation',ch10,&
169 0 : 'of the frequency-dependent dielectric tensor.',ch10,&
170 0 : 'Action : likely there is no action to be taken, although you,',ch10,&
171 0 : 'could try to raise your convergence parameters (ecut and k-points).',ch10,&
172 0 : 'For your information, here are the four lowest frequencies :',ch10,&
173 0 : (phfrq(ii),ii=1,4)
174 0 : ABI_ERROR(message)
175 : end if
176 :
177 : ! Compute the relaxed ion dielectric tensor
178 40 : do idir1=1,3
179 130 : do idir2=1,3
180 : ! The electronic contribution to epsilon is added
181 90 : dielt_rlx(idir1,idir2)=epsinf(idir1,idir2)
182 : ! calculation of the phonon contribution (ionic) to epsilon
183 1227 : do imode=4,3*natom
184 : ! Note that the acoustic modes are not included : their
185 : ! oscillator strength should be exactly zero
186 : ! Also, only the real part of oscstr is taken into account:
187 : ! the possible imaginary parts of degenerate modes
188 : ! will cancel.
189 : dielt_rlx(idir1,idir2)=dielt_rlx(idir1,idir2)+&
190 1107 : & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol
191 : ! Mode decomposition of epsilon
192 1197 : if (dieflag==3)then
193 0 : dielt_modedecompo(idir1,idir2,imode)=oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol
194 : endif ! mode decompo of epsilon
195 : !DEBUG
196 : ! if(idir1==1 .and. idir2==2)then
197 : ! write(std_out,'(a,i4,a,3es16.6)')'imode=',imode,' dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)=',&
198 : !& dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)
199 : ! endif
200 : !ENDDEBUG
201 : end do ! imode
202 : end do ! idir2
203 : end do ! idir1
204 :
205 : ! Print the electronic dielectric tensor
206 10 : call ddb_diel_elec(iout,epsinf)
207 :
208 : ! Print the relaxed ion dielectric tensor
209 10 : write(message,'(a,a)') ch10,' Relaxed ion dielectric tensor'
210 10 : call wrtout(std_out,message,'COLL')
211 10 : call wrtout(iout,message,'COLL')
212 :
213 40 : do idir1=1,3
214 30 : write(message,'(3f16.8)')(dielt_rlx(idir1,idir2),idir2=1,3)
215 30 : call wrtout(std_out,message,'COLL')
216 40 : call wrtout(iout,message,'COLL')
217 : end do
218 10 : call wrtout(iout, " ",'COLL')
219 10 : call wrtout(std_out, " ",'COLL')
220 :
221 : ! Mode decompo of epsilon
222 10 : if (dieflag==3) then
223 0 : write(message,'(a,a,a,a)') ch10,' Mode by mode decomposition of the ionic dielectric tensor',&
224 0 : ch10,' (the electronic contribution is not included)'
225 0 : call wrtout(std_out,message,'COLL')
226 0 : call wrtout(iout,message,'COLL')
227 0 : do imode=4,3*natom
228 0 : write(message,'(a,a,i4,a,es14.6,a,a,3f8.3)') ch10,' Mode number ',imode, ' freq = ',phfrq(imode),' Ha', &
229 0 : ' Mode Z* (|x|, |y|, |z|)', (abs(modez(1,idir1,imode)),idir1=1,3)
230 0 : call wrtout(std_out,message,'COLL')
231 0 : call wrtout(iout,message,'COLL')
232 0 : do idir1=1,3
233 : ! write(message,'(a,a,i4)') ch10,' Mode number 2',imode
234 0 : write(message,'(3f16.8)')(dielt_modedecompo(idir1,idir2,imode),idir2=1,3)
235 0 : call wrtout(std_out,message,'COLL')
236 0 : call wrtout(iout,message,'COLL')
237 : end do
238 : end do
239 : endif ! mode decompo of epsilon
240 :
241 : ! write the relaxed ion dielectric tensor to the netcdf
242 10 : if (ncid /= nctk_noid) then
243 : ncerr = nctk_def_arrays(ncid, [nctkarr_t("emacro_cart_rlx", "dp", &
244 20 : "number_of_cartesian_directions, number_of_cartesian_directions")],defmode=.True.)
245 10 : NCF_CHECK(ncerr)
246 10 : NCF_CHECK(nctk_set_datamode(ncid))
247 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart_rlx"), dielt_rlx))
248 : end if
249 :
250 : ! Frequency-dependent dielectric tensor case
251 10 : if(dieflag==1) then
252 8 : write(message,'(4a)') ch10,ch10,' Calculate the freq. dep. dielectric constant',ch10
253 8 : call wrtout(std_out,message,'COLL')
254 8 : write(message,'(3a)') ch10,' Frequency dependent dielectric constant:',ch10
255 8 : call wrtout(iout,message,'COLL')
256 :
257 : ! Check the possibility of asking the frequency-dependent
258 : ! dielectric tensor (there should be more than one atom in the unit cell)
259 : ! EB: this check is not in any automatic test
260 8 : if(natom==1)then
261 : write(message, '(6a)' )&
262 0 : ' ddb_diel : WARNING -',ch10,&
263 0 : ' When there is only one atom in the unit cell',ch10,&
264 0 : ' cell, the dielectric tensor is frequency-independent.'!,&
265 : ! ' Consequently, dieflag has been reset to 2 . '
266 0 : call wrtout(std_out,message,'COLL')
267 0 : call wrtout(iout,message,'COLL')
268 : end if
269 :
270 8 : difffr=zero
271 8 : if(nfreq>1)difffr=(anaddb_dtset%frmax-anaddb_dtset%frmin)/(nfreq-1)
272 :
273 2 : if (nfreq>10 .and. my_rank == master) then
274 : write(iout, '(a,a,a,a,a,a,a,a)' )&
275 2 : ' ddb_diel : the number of frequencies is larger',&
276 2 : ' than 10 => I will consider only',ch10,&
277 2 : ' the three principal directions, assume that the tensor',ch10,&
278 2 : ' is diagonalized, and give dielectric constant and ',ch10,&
279 4 : ' reflectivities.'
280 : write(iout, '(a,a)' )&
281 2 : ' Frequency(Hartree) Dielectric constant ',&
282 4 : ' Reflectivity '
283 : write(iout, '(a,a)' )&
284 2 : ' x y z',&
285 4 : ' x y z'
286 : end if
287 :
288 : ! Loop on frequencies
289 214 : do ifreq=1,nfreq
290 206 : afreq=anaddb_dtset%frmin+difffr*(ifreq-1)
291 824 : do idir1=1,3
292 2678 : do idir2=1,3
293 1854 : frdiel(idir1,idir2,ifreq)=epsinf(idir1,idir2)
294 8817 : do imode=4,3*natom
295 : ! Note that the acoustic modes are not included : their
296 : ! oscillator strength should be exactly zero
297 : ! Also, only the real part of oscstr is taken into account:
298 : ! the possible imaginary parts of degenerate modes
299 : ! will cancel.
300 : frdiel(idir1,idir2,ifreq)=frdiel(idir1,idir2,ifreq)+&
301 8199 : & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2-afreq**2)*four_pi/ucvol
302 : end do
303 : end do
304 : end do
305 :
306 : ! Write all this information (actually, there should be a choice of units for the frequencies ...
307 214 : if (nfreq>10) then
308 800 : do idir1=1,3
309 800 : if(frdiel(idir1,idir1,ifreq)<=zero)then
310 24 : refl(idir1)=one
311 : else
312 : ! See Gervais and Piriou PRB11,3944(1975) [[cite:Gervais1975]].
313 576 : refl(idir1)=( (sqrt(frdiel(idir1,idir1,ifreq)) -one) /(sqrt(frdiel(idir1,idir1,ifreq)) +one) )**2
314 : end if
315 : end do
316 200 : if (my_rank == master) then
317 800 : write(iout, '(7es12.4)' )afreq,(frdiel(idir1,idir1,ifreq),idir1=1,3),(refl(idir1),idir1=1,3)
318 : end if
319 :
320 : else
321 6 : if (my_rank == master) then
322 6 : write(iout, '(a,es12.4,a)' )' Full dielectric tensor at frequency',afreq,' Hartree'
323 24 : do idir1=1,3
324 24 : write(iout, '(3f16.8)' ) (frdiel(idir1,idir2,ifreq),idir2=1,3)
325 : end do
326 6 : write(iout, '(a)' )' '
327 : end if
328 : end if
329 :
330 : end do ! End of the loop on frequencies
331 : end if ! End the condition on frequency-dependent dielectric tensor (dieflag=1)
332 : end if ! dieflag/=2 .and. nph2l==0
333 : end if
334 :
335 : !Calculation of the Lyddane-Sachs-Teller value of the dielectric constant at zero frequency
336 26 : if(nph2l/=0 .and. dieflag/=2) then
337 :
338 : ! Prepare the output
339 9 : write(message, '(a,a,a,a)' ) ch10,&
340 9 : ' Generalized Lyddane-Sachs-Teller relation at zero frequency :',ch10,&
341 18 : ' Direction Dielectric constant'
342 9 : call wrtout(std_out,message,'COLL')
343 9 : call wrtout(iout,message,'COLL')
344 :
345 : ! Examine every wavevector in the phonon list
346 23 : do iphl2=1,anaddb_dtset%nph2l
347 :
348 56 : qphon(1:3)=anaddb_dtset%qph2l(1:3,iphl2)
349 :
350 23 : if(abs(qphon(1))>DDB_QTOL .or. abs(qphon(2))>DDB_QTOL .or. abs(qphon(3))>DDB_QTOL)then
351 14 : q2=qphon(1)**2+qphon(2)**2+qphon(3)**2
352 : eps=qphon(1)**2*epsinf(1,1)+qphon(2)**2*epsinf(2,2)+&
353 : & qphon(3)**2*epsinf(3,3)+ 2* ( qphon(1)*qphon(2)*epsinf(1,2)+&
354 14 : & qphon(1)*qphon(3)*epsinf(1,3)+qphon(2)*qphon(3)*epsinf(2,3))
355 14 : eps=eps/q2*exp(lst(iphl2)-lst(anaddb_dtset%nph2l+1))
356 14 : if (my_rank == master) then
357 14 : write(iout, '(3f10.5,f16.8)' )qphon,eps
358 14 : write(std_out,'(3f10.5,f16.8)' )qphon,eps
359 : end if
360 : end if
361 : end do
362 : end if ! End of the condition of nph2l does not vanish for Lyddane-Sachs-Teller
363 :
364 :
365 26 : ABI_FREE(frdiel)
366 26 : ABI_FREE(modez)
367 26 : ABI_FREE(oscstr)
368 26 : ABI_FREE(dielt_modedecompo)
369 :
370 26 : end subroutine ddb_diel
371 : !!***
372 :
373 :
374 :
375 : !!****f* ABINIT/ddb_diel_elec
376 : !!
377 : !! NAME
378 : !! ddb_diel_elec
379 : !!
380 : !! FUNCTION
381 : !! Print the electronic dielectric constant (clamped ions)
382 : !!
383 : !! INPUTS
384 : !! iout=unit number for outputs
385 : !! epsinf(3,3)= epsilon^infty = electronic contribution to the
386 : !! dielectric tensor
387 : !!
388 : !! OUTPUT
389 : !!
390 : !! SOURCE
391 :
392 17 : subroutine ddb_diel_elec(iout,epsinf)
393 :
394 : !Arguments -------------------------------
395 : !scalars
396 : integer,intent(in) :: iout
397 : !arrays
398 : real(dp),intent(in) :: epsinf(3,3)
399 :
400 : !Local variables -------------------------
401 : !scalars
402 : integer :: idir1,idir2
403 : character(len=500) :: message
404 : !arrays
405 :
406 17 : write(message, '(a,a)' ) ch10,' Electronic dielectric tensor'
407 17 : call wrtout(std_out,message,'COLL')
408 17 : call wrtout(iout,message,'COLL')
409 :
410 : !Compute the electronic contribution to the dielectric tensor
411 : !Needs only the perturbations with E-field from the DDB
412 68 : do idir1=1,3
413 : ! do idir2=1,3
414 : ! epsinf(idir1,idir2)=d2cart(1,idir1,natom+2,idir2,natom+2)
415 : ! end do
416 51 : write(message, '(3f16.8)' )(epsinf(idir1,idir2),idir2=1,3)
417 51 : call wrtout(std_out,message,'COLL')
418 68 : call wrtout(iout,message,'COLL')
419 : end do
420 17 : call wrtout(iout, " ",'COLL')
421 17 : call wrtout(std_out, " ",'COLL')
422 :
423 17 : end subroutine ddb_diel_elec
424 : !!***
425 :
426 :
427 :
428 : !!****f* ABINIT/ddb_oscstr
429 : !!
430 : !! NAME
431 : !! ddb_oscstr
432 : !!
433 : !! FUNCTION
434 : !! Compute the oscillator strength and the mode effective charge
435 : !!
436 : !! INPUTS
437 : !! displ(2,3*natom,3*natom)=
438 : !! the displacements of atoms in cartesian coordinates.
439 : !! The first index means either the real or the imaginary part,
440 : !! The second index runs on the direction and the atoms displaced
441 : !! The third index runs on the modes.
442 : !! d2cart(2,3,mpert,3,mpert)=
443 : !! dynamical matrix, effective charges, dielectric tensor,....
444 : !! all in cartesian coordinates
445 : !! iout=unit number for outputs
446 : !! mpert =maximum number of ipert
447 : !! natom=number of atoms in unit cell
448 : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
449 : !! matrix eigenvalues, except if these are negative, and in this
450 : !! case, give minus the square root of the absolute value
451 : !! of the matrix eigenvalues). Hartree units.
452 : !!
453 : !! OUTPUT
454 : !! fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes,
455 : !! for different direction of the electric field.
456 : !! modez(2,3,3*natom)=mode effective charges for the different eigenmodes,
457 : !! for different directions of the electric field, following
458 : !! the definition Eq.(53) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
459 : !! oscstr(2,3,3,3*natom)=oscillator strengths, following
460 : !! the definition Eq.(54) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
461 : !!
462 : !! SOURCE
463 :
464 10 : subroutine ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank)
465 :
466 : !Arguments -------------------------------
467 : !scalars
468 : integer,intent(in) :: iout,mpert,natom,ncid,my_rank
469 : !arrays
470 : real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
471 : real(dp),intent(in) :: phfrq(3*natom)
472 : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
473 : real(dp),intent(out) :: fact_oscstr(2,3,3*natom),oscstr(2,3,3,3*natom),modez(2,3,3*natom)
474 :
475 : !Local variables -------------------------
476 : !scalars
477 : integer,parameter :: master=0
478 : integer :: i1,idir1,idir2,imode,ipert1
479 : integer :: ncerr
480 : real(dp) :: usquare
481 : logical :: t_degenerate
482 : !arrays
483 10 : character(len=1),allocatable :: metacharacter(:)
484 :
485 : ! *********************************************************************
486 :
487 : ! Get the factors of the oscillator strength, and the mode effective charge for each mode
488 163 : do imode=1,3*natom
489 : usquare=zero
490 3438 : do i1=1,3*natom
491 3438 : usquare=usquare+displ(1,i1,imode)*displ(1,i1,imode)+displ(2,i1,imode)*displ(2,i1,imode)
492 : end do
493 622 : do idir2=1,3
494 1377 : fact_oscstr(:,idir2,imode)=zero
495 1377 : modez(:,idir2,imode)=zero
496 1989 : do idir1=1,3
497 11691 : do ipert1=1,natom
498 9855 : i1=idir1+(ipert1-1)*3
499 : fact_oscstr(:,idir2,imode)=fact_oscstr(:,idir2,imode)+&
500 29565 : & displ(:,i1,imode)*d2cart(1,idir1,ipert1,idir2,natom+2)
501 : modez(:,idir2,imode)=modez(:,idir2,imode)+&
502 : & displ(:,i1,imode)*&
503 30942 : & d2cart(1,idir1,ipert1,idir2,natom+2)/sqrt(usquare)
504 : end do
505 : end do
506 : end do
507 : end do
508 :
509 : ! Examine the degeneracy of each mode. The portability of the echo of the mode effective
510 : ! charges and oscillator strengths is very hard to guarantee. On the contrary,
511 : ! the scalar reductions of these quantities are OK.
512 40 : ABI_MALLOC(metacharacter,(3*natom))
513 163 : do imode=1,3*natom
514 : ! The degenerate modes are not portable
515 153 : t_degenerate=.false.
516 153 : if(imode>1)then
517 143 : if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
518 : end if
519 153 : if(imode<3*natom)then
520 143 : if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
521 : end if
522 153 : metacharacter(imode)=';'
523 163 : if(t_degenerate)metacharacter(imode)='-'
524 : end do
525 :
526 10 : if (my_rank == master) then
527 : ! Write the mode effective charge for each mode
528 10 : write(iout, '(a)' )' '
529 10 : write(iout, '(a)' )' Mode effective charges '
530 10 : write(iout, '(a)' )' Mode number. x y z length '
531 163 : do imode=1,3*natom
532 153 : write(iout, '(a,i6,a,4f13.6)' )metacharacter(imode),imode,' ',(modez(1,idir1,imode),idir1=1,3),&
533 316 : & (sqrt(modez(1,1,imode)**2+modez(1,2,imode)**2+modez(1,3,imode)**2))
534 : end do
535 : end if ! master
536 :
537 : ! Get the oscillator strengths
538 163 : do imode=1,3*natom
539 622 : do idir1=1,3
540 1989 : do idir2=1,3
541 : oscstr(1,idir1,idir2,imode)= &
542 : & fact_oscstr(1,idir1,imode)*fact_oscstr(1,idir2,imode) +&
543 1377 : & fact_oscstr(2,idir1,imode)*fact_oscstr(2,idir2,imode)
544 1377 : if(abs(oscstr(1,idir1,idir2,imode))<tol14)oscstr(1,idir1,idir2,imode)=zero
545 :
546 : !DEBUG
547 : ! if(idir1==1 .and. idir2==2)then
548 : ! write(std_out,'(a,i4,a,5es16.6)')'imode=',imode,&
549 : !& ' oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)=',&
550 : !& oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)
551 : ! endif
552 : !ENDDEBUG
553 :
554 : oscstr(2,idir1,idir2,imode)= &
555 : & fact_oscstr(1,idir1,imode)*fact_oscstr(2,idir2,imode) -&
556 1377 : & fact_oscstr(2,idir1,imode)*fact_oscstr(1,idir2,imode)
557 1836 : if(abs(oscstr(2,idir1,idir2,imode))<tol14)oscstr(2,idir1,idir2,imode)=zero
558 : end do
559 : end do
560 : end do
561 :
562 10 : if (my_rank == master) then
563 : ! Write the oscillator strength for each mode
564 10 : write(iout, '(a)' )' '
565 10 : write(iout, '(a)' )' Oscillator strengths (in a.u. ; 1 a.u.=253.2638413 m3/s2). Set to zero if abs()<tol14.'
566 10 : write(iout, '(a)' )' Mode number. xx yy zz xy xz yz trace'
567 163 : do imode=1,3*natom
568 : write(iout, '(a,i4,a,7es12.4)' )&
569 765 : & metacharacter(imode),imode,' Real ',(oscstr(1,idir1,idir1,imode),idir1=1,3),&
570 153 : & oscstr(1,1,2,imode), oscstr(1,1,3,imode),oscstr(1,2,3,imode),&
571 918 : & ((oscstr(1,1,1,imode)+oscstr(1,2,2,imode)+oscstr(1,3,3,imode)))
572 : write(iout, '(a,a,6es12.4)' )&
573 765 : & metacharacter(imode),' Imag ',(oscstr(2,idir1,idir1,imode),idir1=1,3),&
574 928 : & oscstr(2,1,2,imode),oscstr(2,1,3,imode),oscstr(2,2,3,imode)
575 : end do
576 :
577 : ! write the oscillator strength to the netcdf
578 10 : if (ncid /= nctk_noid) then
579 : ncerr = nctk_def_arrays(ncid, [nctkarr_t("oscillator_strength", "dp", &
580 : "complex, number_of_cartesian_directions, number_of_cartesian_directions, number_of_phonon_modes")], &
581 20 : defmode=.True.)
582 10 : NCF_CHECK(ncerr)
583 10 : NCF_CHECK(nctk_set_datamode(ncid))
584 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "oscillator_strength"), oscstr))
585 : end if
586 : end if
587 :
588 10 : ABI_FREE(metacharacter)
589 :
590 10 : end subroutine ddb_oscstr
591 : !!***
592 :
593 :
594 :
595 : !!****f* ABINIT/alignph
596 : !!
597 : !! NAME
598 : !! alignph
599 : !!
600 : !! FUNCTION
601 : !! Construct linear combinations of the phonon eigendisplacements
602 : !! of degenerate modes in order to align the mode effective charges
603 : !! along the axes of the cartesian frame.
604 : !!
605 : !! INPUTS
606 : !! amu(ntypat)=mass of the atoms (atomic mass unit)
607 : !! displ(2,3*natom,3*natom)=
608 : !! the displacements of atoms in cartesian coordinates.
609 : !! The first index means either the real or the imaginary part,
610 : !! The second index runs on the direction and the atoms displaced
611 : !! The third index runs on the modes.
612 : !! d2cart(2,3,mpert,3,mpert)=
613 : !! dynamical matrix, effective charges, dielectric tensor,....
614 : !! all in cartesian coordinates
615 : !! mpert =maximum number of ipert
616 : !! natom=number of atoms in unit cell
617 : !! ntypat=number of types of atoms
618 : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
619 : !! matrix eigenvalues, except if these are negative, and in this
620 : !! case, give minus the square root of the absolute value
621 : !! of the matrix eigenvalues). Hartree units.
622 : !! typat(natom)=integer label of each type of atom (1,2,...)
623 : !!
624 : !! OUTPUT
625 : !! displ(2,3*natom,3*natom)=
626 : !! the displacements of atoms in cartesian coordinates.
627 : !! The eigendisplacements of degenerate modes have been aligned along
628 : !! the cartesian axes.
629 : !!
630 : !! SOURCE
631 :
632 33 : subroutine alignph(amu,displ,d2cart,mpert,natom,ntypat,phfrq,typat, &
633 : silent) !optional
634 :
635 : !Arguments -------------------------------
636 : !scalars
637 : integer,intent(in) :: mpert,natom,ntypat
638 : integer,optional,intent(in) :: silent
639 : !arrays
640 : integer,intent(in) :: typat(natom)
641 : real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),phfrq(3*natom)
642 : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
643 :
644 : !Local variables -------------------------
645 : !scalars
646 : integer,parameter :: master=0
647 : integer :: i1,idir1,idir2,ii,imode,imodex,imodey,imodez,ipert1,silent_
648 : real(dp) :: theta
649 : !arrays
650 33 : integer,allocatable :: deg(:)
651 : real(dp) :: zvec(3,3),zvect(3,3)
652 33 : real(dp),allocatable :: modez(:,:,:),modezabs(:),oscstr(:,:,:),vec(:,:),vect(:,:)
653 :
654 : ! *********************************************************************
655 :
656 : !Reduce verbosity
657 33 : silent_=0
658 33 : if (present(silent)) silent_=silent
659 :
660 : !Get the oscillator strength and mode effective charge for each mode
661 132 : ABI_MALLOC(oscstr,(2,3,3*natom))
662 66 : ABI_MALLOC(modez,(2,3,3*natom))
663 99 : ABI_MALLOC(modezabs,(3*natom))
664 132 : ABI_MALLOC(vec,(3*natom,3))
665 66 : ABI_MALLOC(vect,(3*natom,3))
666 99 : ABI_MALLOC(deg,(3*natom))
667 :
668 33 : if (silent_/=1) then
669 3 : write(std_out,'(a,a)')ch10,' alignph : before modifying the eigenvectors, mode number and mode effective charges :'
670 : end if
671 591 : do imode=1,3*natom
672 558 : modezabs(imode)=zero
673 1674 : do ii=1,2
674 5022 : do idir2=1,3
675 3348 : oscstr(ii,idir2,imode)=zero
676 3348 : modez(ii,idir2,imode)=zero
677 13392 : do idir1=1,3
678 72360 : do ipert1=1,natom
679 58968 : i1=idir1+(ipert1-1)*3
680 : oscstr(ii,idir2,imode)=oscstr(ii,idir2,imode)+&
681 : & displ(ii,i1,imode)*&
682 58968 : & d2cart(1,idir1,ipert1,idir2,natom+2)
683 : modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
684 : & displ(ii,i1,imode)*&
685 : & d2cart(1,idir1,ipert1,idir2,natom+2)*&
686 69012 : & sqrt(amu(typat(ipert1))*amu_emass)
687 : end do
688 : end do
689 4464 : if(abs(modez(ii,idir2,imode))>modezabs(imode))modezabs(imode)=abs(modez(ii,idir2,imode))
690 : end do
691 : end do
692 591 : if (silent_/=1) write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
693 : end do
694 :
695 : !Find degenerate modes with non-zero mode effective charge
696 33 : imode = 0
697 482 : do while (imode < 3*natom)
698 449 : imode = imode + 1
699 482 : if (imode == 3*natom) then
700 10 : deg(imode) = 1
701 439 : else if (abs(phfrq(imode) - phfrq(imode+1)) > tol6 .or. modezabs(imode)<tol8 .or. modezabs(imode+1)<tol8) then
702 : ! Differ by phonon frequency or zero mode effective charge
703 336 : deg(imode) = 1
704 : else
705 103 : deg(imode) = 2
706 103 : if (imode < 3*natom - 1) then
707 83 : if (abs(phfrq(imode) - phfrq(imode+2)) < tol6 .and. modezabs(imode+2)>tol8) then
708 6 : deg(imode) = 3
709 6 : imode = imode + 1
710 : end if
711 : end if
712 103 : imode = imode + 1
713 : end if
714 : end do
715 :
716 :
717 : !In case of a degenerate mode, with non-zero mode effective charge, align the mode effective charge vector along
718 : !the axes of the cartesian frame
719 33 : imode = 1
720 482 : do while (imode <= 3*natom)
721 :
722 449 : if (silent_/=1) then
723 6 : write(std_out,'(a,a,i4,a,i2)')ch10,' Mode number ',imode,' has degeneracy ',deg(imode)
724 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of this mode =',modez(1,:,imode)
725 : end if
726 :
727 449 : if (deg(imode) == 2) then
728 :
729 : ! Optimize on the x direction
730 97 : if (silent_/=1) then
731 0 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
732 : end if
733 97 : if (abs(modez(1,1,imode)) > tol8) then
734 95 : theta = atan(-modez(1,1,imode+1)/modez(1,1,imode))
735 1805 : vec(:,1) = displ(1,:,imode)
736 1805 : vec(:,2) = displ(1,:,imode+1)
737 1805 : displ(1,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2)
738 1805 : displ(1,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2)
739 :
740 : ! MR: Rotate also the imaginary part
741 1805 : vec(:,1) = displ(2,:,imode)
742 1805 : vec(:,2) = displ(2,:,imode+1)
743 1805 : displ(2,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2)
744 1805 : displ(2,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2)
745 :
746 : end if
747 :
748 352 : else if (deg(imode) == 3) then
749 :
750 6 : if (silent_/=1) then
751 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
752 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next-next mode =',modez(1,:,imode+2)
753 : end if
754 :
755 : ! Before mixing them, select the mode-effective charge vectors as being predominently "x", "y" or "z" type.
756 0 : if(abs(modez(1,1,imode))>abs(modez(1,2,imode))-tol12 .and. &
757 6 : & abs(modez(1,1,imode))>abs(modez(1,3,imode))-tol12) then
758 0 : imodex=imode
759 0 : if(abs(modez(1,2,imode+1))>abs(modez(1,3,imode+1))-tol12)then
760 0 : imodey=imode+1 ; imodez=imode+2
761 : else
762 0 : imodez=imode+1 ; imodey=imode+2
763 : end if
764 6 : else if(abs(modez(1,2,imode))>abs(modez(1,1,imode))-tol12 .and. &
765 6 : & abs(modez(1,2,imode))>abs(modez(1,3,imode))-tol12) then
766 0 : imodey=imode
767 0 : if(abs(modez(1,1,imode+1))>abs(modez(1,3,imode+1))-tol12)then
768 0 : imodex=imode+1 ; imodez=imode+2
769 : else
770 0 : imodez=imode+1 ; imodex=imode+2
771 : end if
772 : else
773 6 : imodez=imode
774 6 : if(abs(modez(1,1,imode+1))>abs(modez(1,2,imode+1))-tol12)then
775 0 : imodex=imode+1 ; imodey=imode+2
776 : else
777 6 : imodey=imode+1 ; imodex=imode+2
778 : end if
779 : end if
780 42 : vec(:,1)=displ(1,:,imodex)
781 42 : vec(:,2)=displ(1,:,imodey)
782 42 : vec(:,3)=displ(1,:,imodez)
783 24 : zvec(:,1)=modez(1,:,imodex)
784 24 : zvec(:,2)=modez(1,:,imodey)
785 24 : zvec(:,3)=modez(1,:,imodez)
786 :
787 :
788 : ! Optimize along x : does the first vector has a component along x ?
789 6 : if (abs(zvec(1,1)) > tol8) then
790 : ! Optimize on the (1,2) pair of modes along x
791 6 : theta = atan(-zvec(1,2)/zvec(1,1))
792 6 : zvect(:,:)=zvec(:,:)
793 24 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,2)
794 24 : zvec(:,2) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,2)
795 132 : vect(:,:)=vec(:,:)
796 42 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,2)
797 42 : vec(:,2) = sin(theta)*vect(:,1) + cos(theta)*vect(:,2)
798 : ! Optimize on the (1,3) pair of modes along x
799 6 : theta = atan(-zvec(1,3)/zvec(1,1))
800 6 : zvect(:,:)=zvec(:,:)
801 24 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
802 24 : zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
803 132 : vect(:,:)=vec(:,:)
804 42 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
805 42 : vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
806 6 : if (abs(zvec(2,2)) > tol8) then
807 : ! Optimize on the (2,3) pair of modes along y
808 6 : theta = atan(-zvec(2,3)/zvec(2,2))
809 6 : zvect(:,:)=zvec(:,:)
810 24 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
811 24 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
812 132 : vect(:,:)=vec(:,:)
813 42 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
814 42 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
815 : end if
816 : ! Likely, the remaining is not needed ... because the vectors have been ordered in x, y, and z major component ...
817 : ! Optimize along x : does the second vector has a component along x ?
818 0 : else if(abs(zvec(1,2)) > tol8) then
819 : ! Optimize on the (2,3) pair of modes along x
820 0 : theta = atan(-zvec(1,3)/zvec(1,2))
821 0 : zvect(:,:)=zvec(:,:)
822 0 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
823 0 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
824 0 : vect(:,:)=vec(:,:)
825 0 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
826 0 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
827 : ! Optimize on the (1,3) pair of modes along y
828 0 : if (abs(zvec(2,1)) > tol8) then
829 0 : theta = atan(-zvec(2,3)/zvec(2,1))
830 0 : zvect(:,:)=zvec(:,:)
831 0 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
832 0 : zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
833 0 : vect(:,:)=vec(:,:)
834 0 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
835 0 : vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
836 : end if
837 : ! We are left with the pair of vectors (2,3)
838 0 : else if (abs(zvec(2,2)) > tol8) then
839 : ! Optimize on the (2,3) pair of modes along y
840 0 : theta = atan(-zvec(2,3)/zvec(2,2))
841 0 : zvect(:,:)=zvec(:,:)
842 0 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
843 0 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
844 0 : vect(:,:)=vec(:,:)
845 0 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
846 0 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
847 : end if
848 :
849 42 : displ(1,:,imodex)=vec(:,1)
850 42 : displ(1,:,imodey)=vec(:,2)
851 42 : displ(1,:,imodez)=vec(:,3)
852 :
853 : ! Previous coding, from Marek. Apparently, break the orthogonalization of vectors ...
854 : ! do ii = 1,3
855 : ! coeff(:) = 0._dp
856 : ! if (ii == 1) then
857 : ! jj = 2 ; kk = 3
858 : ! else if (ii == 2) then
859 : ! jj = 1 ; kk = 3
860 : ! else
861 : ! jj = 1 ; kk = 2
862 : ! end if
863 : ! coeff(ii) = 1._dp
864 : ! c1 = modez(1,jj,imode+ii-1)
865 : ! c2 = modez(1,jj,imode+jj-1)
866 : ! c3 = modez(1,jj,imode+kk-1)
867 : ! c4 = modez(1,kk,imode+ii-1)
868 : ! c5 = modez(1,kk,imode+jj-1)
869 : ! c6 = modez(1,kk,imode+kk-1)
870 : ! dtm = c2*c6 - c3*c5
871 : ! if (abs(dtm) > tol8) then
872 : ! coeff(jj) = (c3*c4 - c1*c6)/dtm
873 : ! coeff(kk) = (c1*c5 - c2*c4)/dtm
874 : ! end if
875 : ! mod_ = sqrt(1._dp + coeff(jj)*coeff(jj) + coeff(kk)*coeff(kk))
876 : ! coeff(:) = coeff(:)/mod_
877 : ! displ(1,:,imode+ii-1) = coeff(1)*vec(1,:) + coeff(2)*vec(2,:) + &
878 : !& coeff(3)*vec(3,:)
879 : ! end do
880 :
881 : end if ! if deg mode
882 :
883 449 : imode = imode + deg(imode)
884 :
885 : end do
886 :
887 33 : if (silent_/=1) then
888 3 : write(std_out,'(a,a)')ch10,' alignph : after modifying the eigenvectors, mode number and mode effective charges :'
889 : end if
890 591 : do imode=1,3*natom
891 1674 : do ii=1,2
892 5022 : do idir2=1,3
893 3348 : modez(ii,idir2,imode)=zero
894 14508 : do idir1=1,3
895 72360 : do ipert1=1,natom
896 58968 : i1=idir1+(ipert1-1)*3
897 : modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
898 : & displ(ii,i1,imode)*&
899 : & d2cart(1,idir1,ipert1,idir2,natom+2)*&
900 69012 : & sqrt(amu(typat(ipert1))*amu_emass)
901 : end do
902 : end do
903 : end do
904 : end do
905 591 : if (silent_/=1) write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
906 : end do
907 :
908 33 : ABI_FREE(deg)
909 33 : ABI_FREE(oscstr)
910 33 : ABI_FREE(modez)
911 33 : ABI_FREE(modezabs)
912 33 : ABI_FREE(vec)
913 33 : ABI_FREE(vect)
914 :
915 33 : end subroutine alignph
916 : !!***
917 :
918 : end module m_ddb_diel
919 : !!***
|
