Branch data Line data Source code
1 : : !!****m* ABINIT/m_raman
2 : : !! NAME
3 : : !! m_raman
4 : : !!
5 : : !! FUNCTION
6 : : !! Raman susceptibilities of zone-center phonons and electroo tensor.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (MVeithen)
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_raman
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_nctk
28 : : use netcdf
29 : :
30 : : use m_fstrings, only : sjoin
31 : : use m_matrix, only : matr3inv
32 : :
33 : : implicit none
34 : :
35 : : private
36 : : !!***
37 : :
38 : : public :: ramansus ! Raman susceptibilities of zone-center phonons.
39 : : public :: electrooptic ! Electrooptic tensor and the raman tensors of zone-center phonons.
40 : : public :: defwrite_nonana_raman_terms ! Write raman susceptiblities for q-->0
41 : : public :: defwrite_raman_terms ! Write raman susceptiblities and frequencies for q=0
42 : : !!***
43 : :
44 : : contains
45 : : !!***
46 : :
47 : : !!****f* m_raman/ramansus
48 : : !!
49 : : !! NAME
50 : : !! ramansus
51 : : !!
52 : : !! FUNCTION
53 : : !! Compute the raman susceptibilities of zone-center phonons
54 : : !!
55 : : !! INPUTS
56 : : !! d2cart = second order derivatives of the energy wrt all perturbations
57 : : !! dchide(3,3,3) = non-linear optical coefficients from dtchi
58 : : !! dchidt(natom,3,3,3) = first-order change of the electronic dielectric
59 : : !! tensor induced by an individual atomic displacement
60 : : !! displ = phonon mode atomic displacements
61 : : !! mpert = maximum number of perturbations
62 : : !! natom = number of atoms
63 : : !! phfrq = phonon frequencies
64 : : !! qphnrm=(described below)
65 : : !! ucvol = unit cell volume
66 : : !!
67 : : !! OUTPUT
68 : : !! qphon(3)= to be divided by qphnrm, give the phonon wavevector;
69 : : !! if qphnrm==0.0_dp, then the wavevector is zero (Gamma point)
70 : : !! and qphon gives the direction of
71 : : !! the induced electric field; in the latter case, if qphon is
72 : : !! zero, no non-analytical contribution is included.
73 : : !! rsus
74 : : !!
75 : : !! SOURCE
76 : :
77 : 6 : subroutine ramansus(d2cart,dchide,dchidt,displ,mpert,natom,phfrq,qphon,qphnrm,rsus,ucvol)
78 : :
79 : : !Arguments -----------------------------------
80 : : !scalars
81 : : integer,intent(in) :: mpert,natom
82 : : real(dp),intent(in) :: qphnrm,ucvol
83 : : !arrays
84 : : real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),dchide(3,3,3)
85 : : real(dp),intent(in) :: dchidt(natom,3,3,3),displ(2,3*natom,3*natom)
86 : : real(dp),intent(in) :: phfrq(3*natom)
87 : : real(dp),intent(inout) :: qphon(3)
88 : : real(dp),intent(out) :: rsus(3*natom,3,3)
89 : :
90 : : !Local variables-------------------------------
91 : : !scalars
92 : : integer :: analyt,i1,i1dir,i1pert,i2dir,iatom,idir,imode
93 : : real(dp) :: epsq,fac,g0,g1,g2,qphon2
94 : : logical :: t_degenerate,iwrite
95 : : character(len=500) :: message
96 : : !arrays
97 : : real(dp) :: dijk_q(3,3)
98 : 6 : real(dp),allocatable :: zeff(:,:)
99 : 6 : character(len=1),allocatable :: metacharacter(:)
100 : :
101 : : ! *********************************************************************
102 : :
103 : 6 : iwrite = ab_out > 0
104 : :
105 [ - + + - : 18 : ABI_MALLOC(zeff,(3,natom))
+ - + - -
+ - + ]
106 : :
107 [ + + + + : 402 : rsus(:,:,:) = zero
+ + ]
108 : 6 : epsq = zero
109 [ + + + + ]: 54 : zeff(:,:) = zero
110 : : dijk_q(:,:) = zero
111 : :
112 : : !Determine the analyticity of the matrix.
113 : 6 : analyt=1
114 [ + - ]: 6 : if(abs(qphnrm)<tol8)analyt=0
115 [ + + + + : 6 : if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=1
+ + ]
116 : :
117 : : !In the case the non-analyticity is required :
118 [ + - ]: 4 : if(analyt == 0) then
119 : :
120 : : ! Normalize the limiting direction
121 : 4 : qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2
122 : 4 : qphon(1)=qphon(1)/sqrt(qphon2)
123 : 4 : qphon(2)=qphon(2)/sqrt(qphon2)
124 : 4 : qphon(3)=qphon(3)/sqrt(qphon2)
125 : :
126 : : ! Get the dielectric constant for the limiting direction
127 : : epsq= 0._dp
128 [ + + ]: 16 : do i1dir=1,3
129 [ + + ]: 52 : do i2dir=1,3
130 : 48 : epsq=epsq+qphon(i1dir)*qphon(i2dir)*d2cart(1,i1dir,natom+2,i2dir,natom+2)
131 : : end do
132 : : end do
133 : :
134 : : ! Check if epsq > 0
135 [ - + ]: 4 : if (epsq < tol8) then
136 : 0 : write(message,'(a,es14.6)')' The value of epsq must be > 0 while it is found to be',epsq
137 : 0 : ABI_BUG(message)
138 : : end if
139 : :
140 : : ! Get the effective charges for the limiting direction
141 [ + + ]: 16 : do i1dir=1,3
142 [ + + ]: 40 : do i1pert=1,natom
143 : 24 : zeff(i1dir,i1pert)=zero
144 [ + + ]: 108 : do i2dir=1,3
145 : : zeff(i1dir,i1pert)=zeff(i1dir,i1pert)+qphon(i2dir)*&
146 : 96 : & d2cart(1,i1dir,i1pert,i2dir,natom+2)
147 : : end do
148 : : end do
149 : : end do
150 : :
151 : : ! Get the NLO tensor for the limiting direction !$\sum_{k} d_{ijk} \cdot q_k$
152 : :
153 : 4 : dijk_q(:,:) = zero
154 [ + + ]: 16 : do i1dir = 1, 3
155 [ + + ]: 52 : do i2dir = 1, 3
156 [ + + ]: 156 : do idir = 1, 3
157 : 144 : dijk_q(i1dir,i2dir) = dijk_q(i1dir,i2dir) + dchide(i1dir,i2dir,idir)*qphon(idir)
158 : : end do
159 : : end do
160 : : end do
161 : :
162 : 4 : fac = 16._dp*pi/(ucvol*epsq)
163 [ + + ]: 28 : do imode = 1, 3*natom
164 [ + + ]: 76 : do iatom = 1, natom
165 [ + + ]: 216 : do idir = 1, 3
166 : 144 : i1=idir + (iatom - 1)*3
167 : : rsus(imode,:,:) = rsus(imode,:,:) + &
168 : : & (dchidt(iatom,idir,:,:) - fac*zeff(idir,iatom)*dijk_q(:,:))* &
169 [ + + + + ]: 1920 : & displ(1,i1,imode)
170 : : end do ! disp
171 : : end do ! iatom
172 : : end do ! imode
173 [ + + + + : 268 : rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol)
+ + ]
174 : :
175 : : else
176 [ + + ]: 14 : do imode = 1, 3*natom
177 [ + + ]: 38 : do iatom = 1, natom
178 [ + + ]: 108 : do idir = 1, 3
179 : 72 : i1=idir + (iatom - 1)*3
180 [ + + + + ]: 960 : rsus(imode,:,:) = rsus(imode,:,:) + dchidt(iatom,idir,:,:)*displ(1,i1,imode)
181 : : end do ! disp
182 : : end do ! iatom
183 : : end do ! imode
184 [ + + + + : 134 : rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol)
+ + ]
185 : : end if ! analyt == 0
186 : :
187 [ + + ]: 6 : if (analyt == 0) then
188 [ + - ]: 4 : if (iwrite) then
189 : 4 : write(ab_out,*) ch10
190 : : write(ab_out, '(a,/,a,3f9.5)' )&
191 : 4 : & ' Raman susceptibility of zone-center phonons, with non-analyticity in the',&
192 [ + + ]: 20 : & ' direction (cartesian coordinates)',qphon(1:3)+tol10
193 : : write(ab_out,'(a)')&
194 : 4 : & ' -----------------------------------------------------------------------'
195 : 4 : write(ab_out,*) ch10
196 : : end if
197 : :
198 : : else
199 [ + - ]: 2 : if (iwrite) then
200 : 2 : write(ab_out,*) ch10
201 : 2 : write(ab_out,*)' Raman susceptibilities of transverse zone-center phonon modes'
202 : 2 : write(ab_out,*)' -------------------------------------------------------------'
203 : 2 : write(ab_out,*) ch10
204 : : end if
205 : : end if
206 : :
207 : : !Examine the degeneracy of each mode. The portability of the echo of the Raman susceptibility
208 : : !for each degenerate mode is very hard to guarantee. On the contrary,
209 : : !the scalar reductions of these quantities are OK.
210 [ + - - + : 24 : ABI_MALLOC(metacharacter,(3*natom))
+ - + - -
+ - + ]
211 [ + + ]: 42 : do imode=1,3*natom
212 : : ! The degenerate modes are not portable
213 : 36 : t_degenerate=.false.
214 [ + + ]: 36 : if(imode>1)then
215 [ + + ]: 30 : if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
216 : : end if
217 [ + + ]: 36 : if(imode<3*natom)then
218 [ + + ]: 30 : if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
219 : : end if
220 : 36 : metacharacter(imode)=';'
221 [ + + ]: 42 : if(t_degenerate)metacharacter(imode)='-'
222 : : end do
223 : :
224 [ + + ]: 42 : do imode = 1, 3*natom
225 [ + - ]: 36 : if (iwrite) then
226 : 36 : write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')' Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)'
227 [ + + ]: 144 : do idir = 1,3
228 : 144 : write(ab_out,'(a,4x,3(f16.9,2x))')metacharacter(imode),rsus(imode,idir,:)
229 : : end do
230 : : end if
231 : : ! See R. Caracas and X. Gonze, Thermodynamic Properties of Solids : experiment and modeling, Wiley-VCH,
232 : : ! Ed. S. Chaplot and R. Mittal and N. Choudhury , chap. 8, pp 291-312.
233 : 36 : g0=(rsus(imode,1,1)+rsus(imode,2,2)+rsus(imode,3,3))**2*third
234 : : g1=((rsus(imode,1,2)-rsus(imode,2,1))**2+&
235 : : & (rsus(imode,1,3)-rsus(imode,3,1))**2+&
236 : 36 : & (rsus(imode,2,3)-rsus(imode,3,2))**2)*half
237 : : g2=((rsus(imode,1,2)+rsus(imode,2,1))**2+&
238 : : & (rsus(imode,1,3)+rsus(imode,3,1))**2+&
239 : : & (rsus(imode,2,3)+rsus(imode,3,2))**2)*half +&
240 : : & ((rsus(imode,1,1)-rsus(imode,2,2))**2+&
241 : : & (rsus(imode,2,2)-rsus(imode,3,3))**2+&
242 : 36 : & (rsus(imode,3,3)-rsus(imode,1,1))**2)*third
243 [ + - ]: 42 : if (iwrite) then
244 : 36 : write(ab_out,'(3(a,f16.9))')' Spherical averages : G0=',g0,' G1=',g1,' G2=',g2
245 : 36 : write(ab_out,*)
246 : : end if
247 : : end do
248 : :
249 [ - + ]: 6 : ABI_FREE(metacharacter)
250 [ - + ]: 6 : ABI_FREE(zeff)
251 : :
252 : 6 : end subroutine ramansus
253 : : !!***
254 : :
255 : : !!****f* m_raman/electrooptic
256 : : !! NAME
257 : : !! electrooptic
258 : : !!
259 : : !! FUNCTION
260 : : !! Compute the electrooptic tensor and the raman tensors of zone-center phonons
261 : : !!
262 : : !! INPUTS
263 : : !! dchide(3,3,3) = non-linear optical coefficients
264 : : !! dieflag= dielectric tensor flag. 0=> no dielectric tensor,
265 : : !! 1=> frequency-dependent dielectric tensor,
266 : : !! 2=> only the electronic dielectric tensor.
267 : : !! epsinf=electronic dielectric tensor
268 : : !! fact_oscstr(2,3,3*natom)=factors of the oscillator strengths for the different eigenmodes,
269 : : !! for different direction of the electric field
270 : : !! natom=number of atoms in unit cell
271 : : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
272 : : !! matrix eigenvalues, except if these are negative, and in this
273 : : !! case, give minus the square root of the absolute value
274 : : !! of the matrix eigenvalues). Hartree units.
275 : : !! prtmbm= if equal to 1 write out the mode by mode decomposition of the EO tensor
276 : : !! rsus = Raman susceptibilities
277 : : !! ucvol=unit cell volume
278 : : !!
279 : : !! OUTPUT
280 : : !! (to be completed ?)
281 : : !!
282 : : !! NOTES
283 : : !! 1. The phonon frequencies phfrq should correspond to the
284 : : !! wavevector at Gamma, without any non-analyticities.
285 : : !! 2. Should clean for no imaginary part ...
286 : : !! This routine should be used only by one processor.
287 : : !! 3. frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
288 : : !! mode effective charges for the different eigenmodes,
289 : : !! for different direction of the electric field
290 : : !!
291 : : !! SOURCE
292 : :
293 : 2 : subroutine electrooptic(dchide,dieflag,epsinf,fact_oscstr,natom,phfrq,prtmbm,rsus,ucvol)
294 : :
295 : : !Arguments -------------------------------
296 : : !scalars
297 : : integer,intent(in) :: dieflag,natom,prtmbm
298 : : real(dp),intent(in) :: ucvol
299 : : !arrays
300 : : real(dp),intent(in) :: dchide(3,3,3),epsinf(3,3),fact_oscstr(2,3,3*natom)
301 : : real(dp),intent(in) :: phfrq(3*natom),rsus(3*natom,3,3)
302 : :
303 : : !Local variables -------------------------
304 : : !scalars
305 : : integer :: flag,i1,i2,ii,imode,jj,kk
306 : : real(dp) :: dtm,fac
307 : : logical :: iwrite
308 : : character(len=500) :: message
309 : : !arrays
310 : : integer :: voigtindex(6,2)
311 : : real(dp) :: eta(3,3),rvoigt(6,3),work(3,3,3)
312 : 2 : real(dp),allocatable :: rijk(:,:,:,:),rijk_tot(:,:,:)
313 : :
314 : : ! *********************************************************************
315 : :
316 : : !rijk(1:3*natom,:,:,:) = mode by mode decomposition of the electrooptic tensor
317 : : !rijk(3*natom+1,:,:,:) = electronic contribution
318 : 2 : iwrite = ab_out > 0
319 : :
320 : 2 : voigtindex(1,1) = 1 ; voigtindex(1,2) = 1
321 : 2 : voigtindex(2,1) = 2 ; voigtindex(2,2) = 2
322 : 2 : voigtindex(3,1) = 3 ; voigtindex(3,2) = 3
323 : 2 : voigtindex(4,1) = 2 ; voigtindex(4,2) = 3
324 : 2 : voigtindex(5,1) = 1 ; voigtindex(5,2) = 3
325 : 2 : voigtindex(6,1) = 1 ; voigtindex(6,2) = 2
326 : :
327 [ + - - + : 12 : ABI_MALLOC(rijk,(3*natom+1,3,3,3))
+ - + - +
- + - + -
+ - - + -
+ ]
328 [ - + ]: 2 : ABI_MALLOC(rijk_tot,(3,3,3))
329 [ + + + + : 458 : rijk(:,:,:,:) = 0._dp
+ + + + ]
330 [ + + + + : 80 : rijk_tot(:,:,:) = 0._dp
+ + ]
331 : :
332 : :
333 : : !In case there is no mode with truly negative frequency
334 : : !and the electronic dielectric tensor is available
335 : : !compute the electro-optic tensor
336 : :
337 : 2 : flag = 1
338 : :
339 [ - + ]: 2 : if (abs(phfrq(1)) > abs(phfrq(4))) then
340 : 0 : flag = 0
341 : : write(message,'(6a)')&
342 : 0 : & 'The lowest mode appears to be a "true" negative mode,',ch10,&
343 : 0 : & 'and not an acoustic mode. This precludes the computation',ch10,&
344 : 0 : & 'of the EO tensor.',ch10
345 : 0 : ABI_WARNING(message)
346 : : end if
347 : :
348 : : dtm = epsinf(1,1)*epsinf(2,2)*epsinf(3,3) + &
349 : : & epsinf(1,2)*epsinf(2,3)*epsinf(3,1) + &
350 : : & epsinf(1,3)*epsinf(2,1)*epsinf(3,2) - &
351 : : & epsinf(3,1)*epsinf(2,2)*epsinf(1,3) - &
352 : : & epsinf(3,2)*epsinf(2,3)*epsinf(1,1) - &
353 : 2 : & epsinf(3,3)*epsinf(2,1)*epsinf(1,2)
354 : :
355 [ - + ]: 2 : if (abs(dtm) < tol6) then
356 : 0 : flag = 0
357 : : write(message,'(a,a,a,a,a,a,a,a)')&
358 : 0 : & 'The determinant of the electronic dielectric tensor is zero.',ch10,&
359 : 0 : & 'This preludes the computation fo the EO tensor since',ch10,&
360 : 0 : & 'this quantity requires the inverse of epsilon.',ch10,&
361 : 0 : & 'Action : check you database and the value of dieflag in the input file.',ch10
362 : 0 : ABI_WARNING(message)
363 : : end if
364 : :
365 : : !dieflag is required to be one since the EO tensor
366 : : !requires the oscillator strengths
367 : :
368 [ + - + - ]: 2 : if ((flag == 1).and.(dieflag==1)) then
369 : :
370 : : ! Factor to convert atomic units to MKS units
371 : :
372 : 2 : fac = -16._dp*pi*pi*eps0*(Bohr_Ang**2)*1.0d-8/(e_Cb*sqrt(ucvol))
373 : :
374 : : ! Compute inverse of dielectric tensor
375 : : ! needed to convert the nonlinear optical susceptibility tensor
376 : : ! to the electrooptic tensor
377 : :
378 : 2 : call matr3inv(epsinf,eta)
379 : :
380 [ + - ]: 2 : if (iwrite) then
381 : 2 : write(ab_out,*)ch10
382 : 2 : write(ab_out,*)'Output of the EO tensor (pm/V) in Voigt notations'
383 : 2 : write(ab_out,*)'================================================='
384 : 2 : write(ab_out,*)
385 [ + - ]: 2 : if (prtmbm == 1) then
386 : 2 : write(ab_out,*)'Mode by mode decomposition'
387 : 2 : write(ab_out,*)
388 : : end if
389 : : end if
390 : :
391 : : ! Compute the ionic contribution to the EO tensor
392 : :
393 [ + + ]: 8 : do imode = 4, 3*natom
394 : :
395 [ + - + - ]: 6 : if (prtmbm == 1 .and. iwrite) then
396 : 6 : write(ab_out,*)
397 : 6 : write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')'Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)'
398 : : end if
399 : :
400 [ + + ]: 24 : do ii = 1, 3
401 [ + + ]: 78 : do jj = 1, 3
402 [ + + ]: 234 : do kk = 1, 3
403 : 216 : rijk(imode,ii,jj,kk) = rsus(imode,ii,jj)*fact_oscstr(1,kk,imode)/(phfrq(imode)**2)
404 : : end do
405 : : end do
406 : : end do
407 : :
408 : 6 : work(:,:,:) = 0._dp
409 [ + + ]: 24 : do ii = 1,3
410 [ + + ]: 78 : do jj = 1, 3
411 [ + + ]: 234 : do kk = 1, 3
412 : :
413 [ + + ]: 648 : do i1 = 1, 3
414 [ + + ]: 2106 : do i2 = 1, 3
415 : 1944 : work(ii,jj,kk) = work(ii,jj,kk) + eta(ii,i1)*rijk(imode,i1,i2,kk)*eta(i2,jj)
416 : : end do ! i2
417 : : end do ! i1
418 : :
419 : 162 : rijk(imode,ii,jj,kk) = fac*work(ii,jj,kk)
420 : 216 : rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(imode,ii,jj,kk)
421 : : end do
422 : :
423 : : end do
424 : : end do
425 : :
426 [ + - ]: 8 : if (prtmbm == 1) then
427 : 6 : rvoigt(:,:) = 0._dp
428 [ + + ]: 42 : do i1 = 1, 6
429 : 36 : ii = voigtindex(i1,1)
430 : 36 : jj = voigtindex(i1,2)
431 [ + + ]: 144 : do kk = 1, 3
432 : 144 : rvoigt(i1,kk) = (rijk(imode,ii,jj,kk) + rijk(imode,jj,ii,kk))/2._dp
433 : : end do
434 [ + - ]: 42 : if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
435 : : end do
436 : : end if
437 : :
438 : : end do ! imode
439 : :
440 : : ! Compute the electronic contribution to the EO tensor
441 : :
442 [ + - + - ]: 2 : if (prtmbm == 1 .and. iwrite) then
443 : 2 : write(ab_out,*)
444 : 2 : write(ab_out,*)'Electronic contribution to the EO tensor'
445 : : end if
446 : :
447 : 8 : fac = 16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
448 : :
449 [ + + ]: 8 : do ii = 1,3
450 [ + + ]: 26 : do jj = 1, 3
451 [ + + ]: 78 : do kk = 1, 3
452 : :
453 [ + + ]: 216 : do i1 = 1, 3
454 [ + + ]: 702 : do i2 = 1, 3
455 : : rijk(3*natom+1,ii,jj,kk) = rijk(3*natom+1,ii,jj,kk) + &
456 : 648 : & eta(ii,i1)*dchide(i1,i2,kk)*eta(i2,jj)
457 : : end do ! i2
458 : : end do ! i1
459 : :
460 : 54 : rijk(3*natom+1,ii,jj,kk) = -4._dp*rijk(3*natom+1,ii,jj,kk)*fac
461 : 72 : rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(3*natom+1,ii,jj,kk)
462 : :
463 : : end do
464 : :
465 : : end do
466 : : end do
467 : :
468 [ + - ]: 2 : if (prtmbm == 1) then
469 : 2 : rvoigt(:,:) = 0._dp
470 [ + + ]: 14 : do i1 = 1, 6
471 : 12 : ii = voigtindex(i1,1)
472 : 12 : jj = voigtindex(i1,2)
473 [ + + ]: 48 : do kk = 1, 3
474 : 48 : rvoigt(i1,kk) = (rijk(3*natom+1,ii,jj,kk) + rijk(3*natom+1,jj,ii,kk))/2._dp
475 : : end do
476 [ + - ]: 14 : if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
477 : : end do
478 [ + - ]: 2 : if (iwrite) write(ab_out,*)ch10
479 : : end if
480 : :
481 [ + - ]: 2 : if (iwrite) write(ab_out,*)'Total EO tensor (pm/V) in Voigt notations'
482 : 2 : rvoigt(:,:) = 0._dp
483 [ + + ]: 14 : do i1 = 1, 6
484 : 12 : ii = voigtindex(i1,1)
485 : 12 : jj = voigtindex(i1,2)
486 [ + + ]: 48 : do kk = 1, 3
487 : 48 : rvoigt(i1,kk) = (rijk_tot(ii,jj,kk) + rijk_tot(jj,ii,kk))/2._dp
488 : : end do
489 [ + - ]: 14 : if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
490 : : end do
491 : :
492 : : end if ! flag
493 : :
494 [ - + ]: 2 : ABI_FREE(rijk)
495 [ - + ]: 2 : ABI_FREE(rijk_tot)
496 : :
497 : 2 : end subroutine electrooptic
498 : : !!***
499 : :
500 : : !!****f* m_raman/defwrite_nonana_raman_terms
501 : : !! NAME
502 : : !! defwrite_nonana_raman_terms
503 : : !!
504 : : !! FUNCTION
505 : : !! Write the Raman susceptiblities for q-->0 along different directions in the netcdf file.
506 : : !!
507 : : !! INPUTS
508 : : !! ncid=netcdf file id.
509 : : !! iq_dir=Index of the q-point to be written to file.
510 : : !! ndirs=Number of qpoints.
511 : : !! rsus(3*natom,3,3)=List of Raman susceptibilities along the direction corresponding to iq_dir.
512 : : !! natom=Number of atoms
513 : : !!
514 : : !! OUTPUT
515 : : !! Only writing.
516 : : !!
517 : : !! SOURCE
518 : :
519 : 6 : subroutine defwrite_nonana_raman_terms(ncid, iq_dir, ndirs, natom, rsus, mode)
520 : :
521 : : !Arguments ------------------------------------
522 : : !scalars
523 : : integer,intent(in) :: ncid,natom,iq_dir,ndirs
524 : : character(len=*),intent(in) :: mode
525 : : !arrays
526 : : real(dp),intent(in) :: rsus(3*natom,3,3)
527 : :
528 : : !Local variables-------------------------------
529 : : !scalars
530 : : integer :: ncerr, raman_sus_varid
531 : : ! *************************************************************************
532 : :
533 : : ! Fake use of ndirs, to keep it as argument. This should be removed when ndirs will be used.
534 : : if(.false.) ncerr=ndirs
535 : :
536 : 2 : select case (mode)
537 : : case ("define")
538 [ - + ]: 2 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
539 : : ncerr = nctk_def_arrays(ncid, [ nctkarr_t("non_analytical_raman_sus", "dp", &
540 [ + + ]: 4 : "number_of_non_analytical_directions,number_of_phonon_modes,number_of_cartesian_directions,number_of_cartesian_directions")])
541 [ - + ]: 2 : NCF_CHECK(ncerr)
542 : :
543 [ - + ]: 2 : NCF_CHECK(nctk_set_datamode(ncid))
544 : :
545 : : case ("write")
546 [ - + ]: 4 : NCF_CHECK(nf90_inq_varid(ncid, "non_analytical_raman_sus", raman_sus_varid))
547 [ + + + + ]: 36 : ncerr = nf90_put_var(ncid,raman_sus_varid,rsus, start=[iq_dir,1,1,1], count=[1,3*natom,3,3])
548 [ - + ]: 4 : NCF_CHECK(ncerr)
549 : :
550 : : case default
551 [ + + - ]: 6 : ABI_ERROR(sjoin("Wrong value for mode", mode))
552 : : end select
553 : :
554 : 6 : end subroutine defwrite_nonana_raman_terms
555 : : !!***
556 : :
557 : : !!****f* m_raman/defwrite_raman_terms
558 : : !! NAME
559 : : !! defwrite_raman_terms
560 : : !!
561 : : !! FUNCTION
562 : : !! Write the Raman susceptiblities for q=0 and also the phonon frequncies at gamma.
563 : : !!
564 : : !! INPUTS
565 : : !! ncid=netcdf file id.
566 : : !! rsus(3*natom,3,3)=List of Raman susceptibilities.
567 : : !! natom=Number of atoms
568 : : !!
569 : : !! OUTPUT
570 : : !! Only writing.
571 : : !!
572 : : !! SOURCE
573 : :
574 : 2 : subroutine defwrite_raman_terms(ncid, natom, rsus, phfrq)
575 : :
576 : : !Arguments ------------------------------------
577 : : !scalars
578 : : integer,intent(in) :: ncid,natom
579 : : !arrays
580 : : real(dp),intent(in) :: rsus(3*natom,3,3)
581 : : real(dp),intent(in) :: phfrq(3*natom)
582 : :
583 : : !Local variables-------------------------------
584 : : !scalars
585 : : integer :: ncerr, raman_sus_varid, phmodes_varid
586 : : ! *************************************************************************
587 : :
588 [ - + ]: 2 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
589 : : ncerr = nctk_def_arrays(ncid, [ nctkarr_t("raman_sus", "dp", &
590 : : "number_of_phonon_modes,number_of_cartesian_directions,number_of_cartesian_directions"), &
591 [ + + ]: 6 : nctkarr_t("gamma_phonon_modes", "dp", "number_of_phonon_modes")])
592 [ - + ]: 2 : NCF_CHECK(ncerr)
593 : :
594 [ - + ]: 2 : NCF_CHECK(nctk_set_datamode(ncid))
595 : :
596 [ - + ]: 2 : NCF_CHECK(nf90_inq_varid(ncid, "raman_sus", raman_sus_varid))
597 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid,raman_sus_varid,rsus))
598 [ - + ]: 2 : NCF_CHECK(nf90_inq_varid(ncid, "gamma_phonon_modes", phmodes_varid))
599 [ + - + + : 14 : NCF_CHECK(nf90_put_var(ncid,phmodes_varid,phfrq*Ha_eV))
- + - - ]
600 : :
601 : 2 : end subroutine defwrite_raman_terms
602 : : !!***
603 : :
604 : : end module m_raman
605 : : !!***
|
