TABLE OF CONTENTS


ABINIT/m_raman [ Modules ]

[ Top ] [ Modules ]

NAME

  m_raman

FUNCTION

  Raman susceptibilities of zone-center phonons and electroo tensor.

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (MVeithen)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_raman
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_symtk,           only : matr3inv
34 
35  implicit none
36 
37  private

m_raman/electrooptic [ Functions ]

[ Top ] [ m_raman ] [ Functions ]

NAME

 electrooptic

FUNCTION

 Compute the electrooptic tensor and the raman tensors of zone-center phonons

INPUTS

 dchide(3,3,3) = non-linear optical coefficients
 dieflag= dielectric tensor flag. 0=> no dielectric tensor,
  1=> frequency-dependent dielectric tensor,
  2=> only the electronic dielectric tensor.
 epsinf=electronic dielectric tensor
 fact_oscstr(2,3,3*natom)=factors of the oscillator strengths for the different eigenmodes,
  for different direction of the electric field
 natom=number of atoms in unit cell
 phfrq(3*natom)=phonon frequencies (square root of the dynamical
  matrix eigenvalues, except if these are negative, and in this
  case, give minus the square root of the absolute value
  of the matrix eigenvalues). Hartree units.
 prtmbm= if equal to 1 write out the mode by mode decomposition of the EO tensor
 rsus = Raman susceptibilities
 ucvol=unit cell volume

OUTPUT

  (to be completed ?)

NOTES

 1. The phonon frequencies phfrq should correspond to the
 wavevector at Gamma, without any non-analyticities.
 2. Should clean for no imaginary part ...
 This routine should be used only by one processor.
 3. frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
 mode effective charges for the different eigenmodes,
 for different direction of the electric field

PARENTS

      anaddb

CHILDREN

      matr3inv

SOURCE

313 subroutine electrooptic(dchide,dieflag,epsinf,fact_oscstr,natom,phfrq,prtmbm,rsus,ucvol)
314 
315 
316 !This section has been created automatically by the script Abilint (TD).
317 !Do not modify the following lines by hand.
318 #undef ABI_FUNC
319 #define ABI_FUNC 'electrooptic'
320 !End of the abilint section
321 
322  implicit none
323 
324 !Arguments -------------------------------
325 !scalars
326  integer,intent(in) :: dieflag,natom,prtmbm
327  real(dp),intent(in) :: ucvol
328 !arrays
329  real(dp),intent(in) :: dchide(3,3,3),epsinf(3,3),fact_oscstr(2,3,3*natom)
330  real(dp),intent(in) :: phfrq(3*natom),rsus(3*natom,3,3)
331 
332 !Local variables -------------------------
333 !scalars
334  integer :: flag,i1,i2,ii,imode,jj,kk
335  real(dp) :: dtm,fac
336  logical :: iwrite
337  character(len=500) :: message
338 !arrays
339  integer :: voigtindex(6,2)
340  real(dp) :: eta(3,3),rvoigt(6,3),work(3,3,3)
341  real(dp),allocatable :: rijk(:,:,:,:),rijk_tot(:,:,:)
342 
343 ! *********************************************************************
344 
345 !rijk(1:3*natom,:,:,:) = mode by mode decomposition of the electrooptic tensor
346 !rijk(3*natom+1,:,:,:) = electronic contribution
347  iwrite = ab_out > 0
348 
349  voigtindex(1,1) = 1 ; voigtindex(1,2) = 1
350  voigtindex(2,1) = 2 ; voigtindex(2,2) = 2
351  voigtindex(3,1) = 3 ; voigtindex(3,2) = 3
352  voigtindex(4,1) = 2 ; voigtindex(4,2) = 3
353  voigtindex(5,1) = 1 ; voigtindex(5,2) = 3
354  voigtindex(6,1) = 1 ; voigtindex(6,2) = 2
355 
356  ABI_ALLOCATE(rijk,(3*natom+1,3,3,3))
357  ABI_ALLOCATE(rijk_tot,(3,3,3))
358  rijk(:,:,:,:) = 0._dp
359  rijk_tot(:,:,:) = 0._dp
360 
361 
362 !In case there is no mode with truly negative frequency
363 !and the electronic dielectric tensor is available
364 !compute the electro-optic tensor
365 
366  flag = 1
367 
368  if (abs(phfrq(1)) > abs(phfrq(4))) then
369    flag = 0
370    write(message,'(6a)')&
371 &   'The lowest mode appears to be a "true" negative mode,',ch10,&
372 &   'and not an acoustic mode. This precludes the computation',ch10,&
373 &   'of the EO tensor.',ch10
374    MSG_WARNING(message)
375  end if
376 
377  dtm = epsinf(1,1)*epsinf(2,2)*epsinf(3,3) + &
378 & epsinf(1,2)*epsinf(2,3)*epsinf(3,1) + &
379 & epsinf(1,3)*epsinf(2,1)*epsinf(3,2) - &
380 & epsinf(3,1)*epsinf(2,2)*epsinf(1,3) - &
381 & epsinf(3,2)*epsinf(2,3)*epsinf(1,1) - &
382 & epsinf(3,3)*epsinf(2,1)*epsinf(1,2)
383 
384  if (abs(dtm) < tol6) then
385    flag = 0
386    write(message,'(a,a,a,a,a,a,a,a)')&
387 &   'The determinant of the electronic dielectric tensor is zero.',ch10,&
388 &   'This preludes the computation fo the EO tensor since',ch10,&
389 &   'this quantity requires the inverse of epsilon.',ch10,&
390 &   'Action : check you database and the value of dieflag in the input file.',ch10
391    MSG_WARNING(message)
392  end if
393 
394 !dieflag is required to be one since the EO tensor
395 !requires the oscillator strengths
396 
397  if ((flag == 1).and.(dieflag==1)) then
398 
399 !  Factor to convert atomic units to MKS units
400 
401    fac = -16._dp*pi*pi*eps0*(Bohr_Ang**2)*1.0d-8/(e_Cb*sqrt(ucvol))
402 
403 !  Compute inverse of dielectric tensor
404 !  needed to convert the nonlinear optical susceptibility tensor
405 !  to the electrooptic tensor
406 
407    call matr3inv(epsinf,eta)
408 
409    if (iwrite) then
410      write(ab_out,*)ch10
411      write(ab_out,*)'Output of the EO tensor (pm/V) in Voigt notations'
412      write(ab_out,*)'================================================='
413      write(ab_out,*)
414      if (prtmbm == 1) then
415        write(ab_out,*)'Mode by mode decomposition'
416        write(ab_out,*)
417      end if
418    end if
419 
420 !  Compute the ionic contribution to the EO tensor
421 
422    do imode = 4, 3*natom
423 
424      if (prtmbm == 1 .and. iwrite) then
425        write(ab_out,*)
426        write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')'Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)'
427      end if
428 
429      do ii = 1, 3
430        do jj = 1, 3
431          do kk = 1, 3
432            rijk(imode,ii,jj,kk) = rsus(imode,ii,jj)*fact_oscstr(1,kk,imode)/(phfrq(imode)**2)
433          end do
434        end do
435      end do
436 
437      work(:,:,:) = 0._dp
438      do ii = 1,3
439        do jj = 1, 3
440          do kk = 1, 3
441 
442            do i1 = 1, 3
443              do i2 = 1, 3
444                work(ii,jj,kk) = work(ii,jj,kk) + eta(ii,i1)*rijk(imode,i1,i2,kk)*eta(i2,jj)
445              end do  ! i2
446            end do   ! i1
447 
448            rijk(imode,ii,jj,kk) = fac*work(ii,jj,kk)
449            rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(imode,ii,jj,kk)
450          end do
451 
452        end do
453      end do
454 
455      if (prtmbm == 1) then
456        rvoigt(:,:) = 0._dp
457        do i1 = 1, 6
458          ii = voigtindex(i1,1)
459          jj = voigtindex(i1,2)
460          do kk = 1, 3
461            rvoigt(i1,kk) = (rijk(imode,ii,jj,kk) + rijk(imode,jj,ii,kk))/2._dp
462          end do
463          if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
464        end do
465      end if
466 
467    end do     ! imode
468 
469 !  Compute the electronic contribution to the EO tensor
470 
471    if (prtmbm == 1 .and. iwrite) then
472      write(ab_out,*)
473      write(ab_out,*)'Electronic contribution to the EO tensor'
474    end if
475 
476    fac = 16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
477 
478    do ii = 1,3
479      do jj = 1, 3
480        do kk = 1, 3
481 
482          do i1 = 1, 3
483            do i2 = 1, 3
484              rijk(3*natom+1,ii,jj,kk) = rijk(3*natom+1,ii,jj,kk) + &
485 &             eta(ii,i1)*dchide(i1,i2,kk)*eta(i2,jj)
486            end do  ! i2
487          end do   ! i1
488 
489          rijk(3*natom+1,ii,jj,kk) = -4._dp*rijk(3*natom+1,ii,jj,kk)*fac
490          rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(3*natom+1,ii,jj,kk)
491 
492        end do
493 
494      end do
495    end do
496 
497    if (prtmbm == 1) then
498      rvoigt(:,:) = 0._dp
499      do i1 = 1, 6
500        ii = voigtindex(i1,1)
501        jj = voigtindex(i1,2)
502        do kk = 1, 3
503          rvoigt(i1,kk) = (rijk(3*natom+1,ii,jj,kk) + rijk(3*natom+1,jj,ii,kk))/2._dp
504        end do
505        if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
506      end do
507      if (iwrite) write(ab_out,*)ch10
508    end if
509 
510    if (iwrite) write(ab_out,*)'Total EO tensor (pm/V) in Voigt notations'
511    rvoigt(:,:) = 0._dp
512    do i1 = 1, 6
513      ii = voigtindex(i1,1)
514      jj = voigtindex(i1,2)
515      do kk = 1, 3
516        rvoigt(i1,kk) = (rijk_tot(ii,jj,kk) + rijk_tot(jj,ii,kk))/2._dp
517      end do
518      if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:)
519    end do
520 
521  end if  ! flag
522 
523  ABI_DEALLOCATE(rijk)
524  ABI_DEALLOCATE(rijk_tot)
525 
526 end subroutine electrooptic

m_raman/ramansus [ Functions ]

[ Top ] [ m_raman ] [ Functions ]

NAME

 ramansus

FUNCTION

 Compute the raman susceptibilities of zone-center phonons

INPUTS

  d2cart = second order derivatives of the energy wrt all perturbations
  dchide(3,3,3) = non-linear optical coefficients from dtchi
  dchidt(natom,3,3,3) = first-order change of the electronic dielectric
     tensor induced by an individual atomic displacement
  displ = phonon mode atomic displacements
  mpert = maximum number of perturbations
  natom = number of atoms
  phfrq = phonon frequencies
  qphnrm=(described below)
  ucvol = unit cell volume

OUTPUT

  qphon(3)= to be divided by qphnrm, give the phonon wavevector;
     if qphnrm==0.0_dp, then the wavevector is zero (Gamma point)
     and qphon gives the direction of
     the induced electric field; in the latter case, if qphon is
     zero, no non-analytical contribution is included.
  rsus

PARENTS

      anaddb

CHILDREN

SOURCE

 82 subroutine ramansus(d2cart,dchide,dchidt,displ,mpert,natom,phfrq,qphon,qphnrm,rsus,ucvol)
 83 
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'ramansus'
 89 !End of the abilint section
 90 
 91  implicit none
 92 
 93 !Arguments -----------------------------------
 94 !scalars
 95  integer,intent(in) :: mpert,natom
 96  real(dp),intent(in) :: qphnrm,ucvol
 97 !arrays
 98  real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),dchide(3,3,3)
 99  real(dp),intent(in) :: dchidt(natom,3,3,3),displ(2,3*natom,3*natom)
100  real(dp),intent(in) :: phfrq(3*natom)
101  real(dp),intent(inout) :: qphon(3)
102  real(dp),intent(out) :: rsus(3*natom,3,3)
103 
104 !Local variables-------------------------------
105 !scalars
106  integer :: analyt,i1,i1dir,i1pert,i2dir,iatom,idir,imode
107  real(dp) :: epsq,fac,g0,g1,g2,qphon2
108  logical :: t_degenerate,iwrite
109  character(len=500) :: message
110 !arrays
111  real(dp) :: dijk_q(3,3)
112  real(dp),allocatable :: zeff(:,:)
113  character(len=1),allocatable :: metacharacter(:)
114 
115 ! *********************************************************************
116 
117  iwrite = ab_out > 0
118 
119  ABI_ALLOCATE(zeff,(3,natom))
120 
121  rsus(:,:,:) = zero
122  epsq        = zero
123  zeff(:,:)   = zero
124  dijk_q(:,:) = zero
125 
126 !Determine the analyticity of the matrix.
127  analyt=1
128  if(abs(qphnrm)<tol8)analyt=0
129  if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=1
130 
131 !In the case the non-analyticity is required :
132  if(analyt == 0) then
133 
134 !  Normalize the limiting direction
135    qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2
136    qphon(1)=qphon(1)/sqrt(qphon2)
137    qphon(2)=qphon(2)/sqrt(qphon2)
138    qphon(3)=qphon(3)/sqrt(qphon2)
139 
140 !  Get the dielectric constant for the limiting direction
141    epsq= 0._dp
142    do i1dir=1,3
143      do i2dir=1,3
144        epsq=epsq+qphon(i1dir)*qphon(i2dir)*d2cart(1,i1dir,natom+2,i2dir,natom+2)
145      end do
146    end do
147 
148 !  Check if epsq > 0
149    if (epsq < tol8) then
150      write(message,'(a,es14.6)')'  The value of epsq must be > 0 while it is found to be',epsq
151      MSG_BUG(message)
152    end if
153 
154 !  Get the effective charges for the limiting direction
155    do i1dir=1,3
156      do i1pert=1,natom
157        zeff(i1dir,i1pert)=zero
158        do i2dir=1,3
159          zeff(i1dir,i1pert)=zeff(i1dir,i1pert)+qphon(i2dir)*&
160 &         d2cart(1,i1dir,i1pert,i2dir,natom+2)
161        end do
162      end do
163    end do
164 
165 !  Get the NLO tensor for the limiting direction !$\sum_{k} d_{ijk} \cdot q_k$
166 
167    dijk_q(:,:) = zero
168    do i1dir = 1, 3
169      do i2dir = 1, 3
170        do idir = 1, 3
171          dijk_q(i1dir,i2dir) = dijk_q(i1dir,i2dir) + dchide(i1dir,i2dir,idir)*qphon(idir)
172        end do
173      end do
174    end do
175 
176    fac = 16._dp*pi/(ucvol*epsq)
177    do imode = 1, 3*natom
178      do iatom = 1, natom
179        do idir = 1, 3
180          i1=idir + (iatom - 1)*3
181          rsus(imode,:,:) = rsus(imode,:,:) + &
182 &         (dchidt(iatom,idir,:,:) - fac*zeff(idir,iatom)*dijk_q(:,:))* &
183 &         displ(1,i1,imode)
184        end do  ! disp
185      end do  ! iatom
186    end do  ! imode
187    rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol)
188 
189  else
190    do imode = 1, 3*natom
191      do iatom = 1, natom
192        do idir = 1, 3
193          i1=idir + (iatom - 1)*3
194          rsus(imode,:,:) = rsus(imode,:,:) + dchidt(iatom,idir,:,:)*displ(1,i1,imode)
195        end do  ! disp
196      end do  ! iatom
197    end do  ! imode
198    rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol)
199  end if      ! analyt == 0
200 
201  if (analyt == 0) then
202    if (iwrite) then
203      write(ab_out,*) ch10
204      write(ab_out, '(a,/,a,3f9.5)' )&
205 &     ' Raman susceptibility of zone-center phonons, with non-analyticity in the',&
206 &     '  direction (cartesian coordinates)',qphon(1:3)+tol10
207      write(ab_out,'(a)')&
208 &     ' -----------------------------------------------------------------------'
209      write(ab_out,*) ch10
210    end if
211 
212  else
213    if (iwrite) then
214      write(ab_out,*) ch10
215      write(ab_out,*)' Raman susceptibilities of transverse zone-center phonon modes'
216      write(ab_out,*)' -------------------------------------------------------------'
217      write(ab_out,*) ch10
218    end if
219  end if
220 
221 !Examine the degeneracy of each mode. The portability of the echo of the Raman susceptibility
222 !for each degenerate mode is very hard to guarantee. On the contrary,
223 !the scalar reductions of these quantities are OK.
224  ABI_ALLOCATE(metacharacter,(3*natom))
225  do imode=1,3*natom
226 !  The degenerate modes are not portable
227    t_degenerate=.false.
228    if(imode>1)then
229      if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
230    end if
231    if(imode<3*natom)then
232      if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
233    end if
234    metacharacter(imode)=';'
235    if(t_degenerate)metacharacter(imode)='-'
236  end do
237 
238  do imode = 1, 3*natom
239    if (iwrite) then
240      write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')' Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)'
241      do idir = 1,3
242        write(ab_out,'(a,4x,3(f16.9,2x))')metacharacter(imode),rsus(imode,idir,:)
243      end do
244    end if
245 !  See R. Caracas and X. Gonze, Thermodynamic Properties of Solids : experiment and modeling, Wiley-VCH,
246 !  Ed. S. Chaplot and R. Mittal and N. Choudhury , chap. 8, pp 291-312.
247    g0=(rsus(imode,1,1)+rsus(imode,2,2)+rsus(imode,3,3))**2*third
248    g1=((rsus(imode,1,2)-rsus(imode,2,1))**2+&
249 &   (rsus(imode,1,3)-rsus(imode,3,1))**2+&
250 &   (rsus(imode,2,3)-rsus(imode,3,2))**2)*half
251    g2=((rsus(imode,1,2)+rsus(imode,2,1))**2+&
252 &   (rsus(imode,1,3)+rsus(imode,3,1))**2+&
253 &   (rsus(imode,2,3)+rsus(imode,3,2))**2)*half +&
254 &   ((rsus(imode,1,1)-rsus(imode,2,2))**2+&
255 &   (rsus(imode,2,2)-rsus(imode,3,3))**2+&
256 &   (rsus(imode,3,3)-rsus(imode,1,1))**2)*third
257    if (iwrite) then
258      write(ab_out,'(3(a,f16.9))')' Spherical averages : G0=',g0,'    G1=',g1,'    G2=',g2
259      write(ab_out,*)
260    end if
261  end do
262 
263  ABI_DEALLOCATE(metacharacter)
264  ABI_DEALLOCATE(zeff)
265 
266 end subroutine ramansus