TABLE OF CONTENTS
ABINIT/m_raman [ Modules ]
NAME
m_raman
FUNCTION
Raman susceptibilities of zone-center phonons and electroo tensor.
COPYRIGHT
Copyright (C) 1999-2024 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 .
SOURCE
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 28 use m_symtk, only : matr3inv 29 30 implicit none 31 32 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
SOURCE
288 subroutine electrooptic(dchide,dieflag,epsinf,fact_oscstr,natom,phfrq,prtmbm,rsus,ucvol) 289 290 !Arguments ------------------------------- 291 !scalars 292 integer,intent(in) :: dieflag,natom,prtmbm 293 real(dp),intent(in) :: ucvol 294 !arrays 295 real(dp),intent(in) :: dchide(3,3,3),epsinf(3,3),fact_oscstr(2,3,3*natom) 296 real(dp),intent(in) :: phfrq(3*natom),rsus(3*natom,3,3) 297 298 !Local variables ------------------------- 299 !scalars 300 integer :: flag,i1,i2,ii,imode,jj,kk 301 real(dp) :: dtm,fac 302 logical :: iwrite 303 character(len=500) :: message 304 !arrays 305 integer :: voigtindex(6,2) 306 real(dp) :: eta(3,3),rvoigt(6,3),work(3,3,3) 307 real(dp),allocatable :: rijk(:,:,:,:),rijk_tot(:,:,:) 308 309 ! ********************************************************************* 310 311 !rijk(1:3*natom,:,:,:) = mode by mode decomposition of the electrooptic tensor 312 !rijk(3*natom+1,:,:,:) = electronic contribution 313 iwrite = ab_out > 0 314 315 voigtindex(1,1) = 1 ; voigtindex(1,2) = 1 316 voigtindex(2,1) = 2 ; voigtindex(2,2) = 2 317 voigtindex(3,1) = 3 ; voigtindex(3,2) = 3 318 voigtindex(4,1) = 2 ; voigtindex(4,2) = 3 319 voigtindex(5,1) = 1 ; voigtindex(5,2) = 3 320 voigtindex(6,1) = 1 ; voigtindex(6,2) = 2 321 322 ABI_MALLOC(rijk,(3*natom+1,3,3,3)) 323 ABI_MALLOC(rijk_tot,(3,3,3)) 324 rijk(:,:,:,:) = 0._dp 325 rijk_tot(:,:,:) = 0._dp 326 327 328 !In case there is no mode with truly negative frequency 329 !and the electronic dielectric tensor is available 330 !compute the electro-optic tensor 331 332 flag = 1 333 334 if (abs(phfrq(1)) > abs(phfrq(4))) then 335 flag = 0 336 write(message,'(6a)')& 337 & 'The lowest mode appears to be a "true" negative mode,',ch10,& 338 & 'and not an acoustic mode. This precludes the computation',ch10,& 339 & 'of the EO tensor.',ch10 340 ABI_WARNING(message) 341 end if 342 343 dtm = epsinf(1,1)*epsinf(2,2)*epsinf(3,3) + & 344 & epsinf(1,2)*epsinf(2,3)*epsinf(3,1) + & 345 & epsinf(1,3)*epsinf(2,1)*epsinf(3,2) - & 346 & epsinf(3,1)*epsinf(2,2)*epsinf(1,3) - & 347 & epsinf(3,2)*epsinf(2,3)*epsinf(1,1) - & 348 & epsinf(3,3)*epsinf(2,1)*epsinf(1,2) 349 350 if (abs(dtm) < tol6) then 351 flag = 0 352 write(message,'(a,a,a,a,a,a,a,a)')& 353 & 'The determinant of the electronic dielectric tensor is zero.',ch10,& 354 & 'This preludes the computation fo the EO tensor since',ch10,& 355 & 'this quantity requires the inverse of epsilon.',ch10,& 356 & 'Action : check you database and the value of dieflag in the input file.',ch10 357 ABI_WARNING(message) 358 end if 359 360 !dieflag is required to be one since the EO tensor 361 !requires the oscillator strengths 362 363 if ((flag == 1).and.(dieflag==1)) then 364 365 ! Factor to convert atomic units to MKS units 366 367 fac = -16._dp*pi*pi*eps0*(Bohr_Ang**2)*1.0d-8/(e_Cb*sqrt(ucvol)) 368 369 ! Compute inverse of dielectric tensor 370 ! needed to convert the nonlinear optical susceptibility tensor 371 ! to the electrooptic tensor 372 373 call matr3inv(epsinf,eta) 374 375 if (iwrite) then 376 write(ab_out,*)ch10 377 write(ab_out,*)'Output of the EO tensor (pm/V) in Voigt notations' 378 write(ab_out,*)'=================================================' 379 write(ab_out,*) 380 if (prtmbm == 1) then 381 write(ab_out,*)'Mode by mode decomposition' 382 write(ab_out,*) 383 end if 384 end if 385 386 ! Compute the ionic contribution to the EO tensor 387 388 do imode = 4, 3*natom 389 390 if (prtmbm == 1 .and. iwrite) then 391 write(ab_out,*) 392 write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')'Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)' 393 end if 394 395 do ii = 1, 3 396 do jj = 1, 3 397 do kk = 1, 3 398 rijk(imode,ii,jj,kk) = rsus(imode,ii,jj)*fact_oscstr(1,kk,imode)/(phfrq(imode)**2) 399 end do 400 end do 401 end do 402 403 work(:,:,:) = 0._dp 404 do ii = 1,3 405 do jj = 1, 3 406 do kk = 1, 3 407 408 do i1 = 1, 3 409 do i2 = 1, 3 410 work(ii,jj,kk) = work(ii,jj,kk) + eta(ii,i1)*rijk(imode,i1,i2,kk)*eta(i2,jj) 411 end do ! i2 412 end do ! i1 413 414 rijk(imode,ii,jj,kk) = fac*work(ii,jj,kk) 415 rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(imode,ii,jj,kk) 416 end do 417 418 end do 419 end do 420 421 if (prtmbm == 1) then 422 rvoigt(:,:) = 0._dp 423 do i1 = 1, 6 424 ii = voigtindex(i1,1) 425 jj = voigtindex(i1,2) 426 do kk = 1, 3 427 rvoigt(i1,kk) = (rijk(imode,ii,jj,kk) + rijk(imode,jj,ii,kk))/2._dp 428 end do 429 if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:) 430 end do 431 end if 432 433 end do ! imode 434 435 ! Compute the electronic contribution to the EO tensor 436 437 if (prtmbm == 1 .and. iwrite) then 438 write(ab_out,*) 439 write(ab_out,*)'Electronic contribution to the EO tensor' 440 end if 441 442 fac = 16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 443 444 do ii = 1,3 445 do jj = 1, 3 446 do kk = 1, 3 447 448 do i1 = 1, 3 449 do i2 = 1, 3 450 rijk(3*natom+1,ii,jj,kk) = rijk(3*natom+1,ii,jj,kk) + & 451 & eta(ii,i1)*dchide(i1,i2,kk)*eta(i2,jj) 452 end do ! i2 453 end do ! i1 454 455 rijk(3*natom+1,ii,jj,kk) = -4._dp*rijk(3*natom+1,ii,jj,kk)*fac 456 rijk_tot(ii,jj,kk) = rijk_tot(ii,jj,kk) + rijk(3*natom+1,ii,jj,kk) 457 458 end do 459 460 end do 461 end do 462 463 if (prtmbm == 1) then 464 rvoigt(:,:) = 0._dp 465 do i1 = 1, 6 466 ii = voigtindex(i1,1) 467 jj = voigtindex(i1,2) 468 do kk = 1, 3 469 rvoigt(i1,kk) = (rijk(3*natom+1,ii,jj,kk) + rijk(3*natom+1,jj,ii,kk))/2._dp 470 end do 471 if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:) 472 end do 473 if (iwrite) write(ab_out,*)ch10 474 end if 475 476 if (iwrite) write(ab_out,*)'Total EO tensor (pm/V) in Voigt notations' 477 rvoigt(:,:) = 0._dp 478 do i1 = 1, 6 479 ii = voigtindex(i1,1) 480 jj = voigtindex(i1,2) 481 do kk = 1, 3 482 rvoigt(i1,kk) = (rijk_tot(ii,jj,kk) + rijk_tot(jj,ii,kk))/2._dp 483 end do 484 if (iwrite) write(ab_out,'(5x,3(2x,f16.9))')rvoigt(i1,:) 485 end do 486 487 end if ! flag 488 489 ABI_FREE(rijk) 490 ABI_FREE(rijk_tot) 491 492 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
SOURCE
72 subroutine ramansus(d2cart,dchide,dchidt,displ,mpert,natom,phfrq,qphon,qphnrm,rsus,ucvol) 73 74 !Arguments ----------------------------------- 75 !scalars 76 integer,intent(in) :: mpert,natom 77 real(dp),intent(in) :: qphnrm,ucvol 78 !arrays 79 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),dchide(3,3,3) 80 real(dp),intent(in) :: dchidt(natom,3,3,3),displ(2,3*natom,3*natom) 81 real(dp),intent(in) :: phfrq(3*natom) 82 real(dp),intent(inout) :: qphon(3) 83 real(dp),intent(out) :: rsus(3*natom,3,3) 84 85 !Local variables------------------------------- 86 !scalars 87 integer :: analyt,i1,i1dir,i1pert,i2dir,iatom,idir,imode 88 real(dp) :: epsq,fac,g0,g1,g2,qphon2 89 logical :: t_degenerate,iwrite 90 character(len=500) :: message 91 !arrays 92 real(dp) :: dijk_q(3,3) 93 real(dp),allocatable :: zeff(:,:) 94 character(len=1),allocatable :: metacharacter(:) 95 96 ! ********************************************************************* 97 98 iwrite = ab_out > 0 99 100 ABI_MALLOC(zeff,(3,natom)) 101 102 rsus(:,:,:) = zero 103 epsq = zero 104 zeff(:,:) = zero 105 dijk_q(:,:) = zero 106 107 !Determine the analyticity of the matrix. 108 analyt=1 109 if(abs(qphnrm)<tol8)analyt=0 110 if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=1 111 112 !In the case the non-analyticity is required : 113 if(analyt == 0) then 114 115 ! Normalize the limiting direction 116 qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2 117 qphon(1)=qphon(1)/sqrt(qphon2) 118 qphon(2)=qphon(2)/sqrt(qphon2) 119 qphon(3)=qphon(3)/sqrt(qphon2) 120 121 ! Get the dielectric constant for the limiting direction 122 epsq= 0._dp 123 do i1dir=1,3 124 do i2dir=1,3 125 epsq=epsq+qphon(i1dir)*qphon(i2dir)*d2cart(1,i1dir,natom+2,i2dir,natom+2) 126 end do 127 end do 128 129 ! Check if epsq > 0 130 if (epsq < tol8) then 131 write(message,'(a,es14.6)')' The value of epsq must be > 0 while it is found to be',epsq 132 ABI_BUG(message) 133 end if 134 135 ! Get the effective charges for the limiting direction 136 do i1dir=1,3 137 do i1pert=1,natom 138 zeff(i1dir,i1pert)=zero 139 do i2dir=1,3 140 zeff(i1dir,i1pert)=zeff(i1dir,i1pert)+qphon(i2dir)*& 141 & d2cart(1,i1dir,i1pert,i2dir,natom+2) 142 end do 143 end do 144 end do 145 146 ! Get the NLO tensor for the limiting direction !$\sum_{k} d_{ijk} \cdot q_k$ 147 148 dijk_q(:,:) = zero 149 do i1dir = 1, 3 150 do i2dir = 1, 3 151 do idir = 1, 3 152 dijk_q(i1dir,i2dir) = dijk_q(i1dir,i2dir) + dchide(i1dir,i2dir,idir)*qphon(idir) 153 end do 154 end do 155 end do 156 157 fac = 16._dp*pi/(ucvol*epsq) 158 do imode = 1, 3*natom 159 do iatom = 1, natom 160 do idir = 1, 3 161 i1=idir + (iatom - 1)*3 162 rsus(imode,:,:) = rsus(imode,:,:) + & 163 & (dchidt(iatom,idir,:,:) - fac*zeff(idir,iatom)*dijk_q(:,:))* & 164 & displ(1,i1,imode) 165 end do ! disp 166 end do ! iatom 167 end do ! imode 168 rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol) 169 170 else 171 do imode = 1, 3*natom 172 do iatom = 1, natom 173 do idir = 1, 3 174 i1=idir + (iatom - 1)*3 175 rsus(imode,:,:) = rsus(imode,:,:) + dchidt(iatom,idir,:,:)*displ(1,i1,imode) 176 end do ! disp 177 end do ! iatom 178 end do ! imode 179 rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol) 180 end if ! analyt == 0 181 182 if (analyt == 0) then 183 if (iwrite) then 184 write(ab_out,*) ch10 185 write(ab_out, '(a,/,a,3f9.5)' )& 186 & ' Raman susceptibility of zone-center phonons, with non-analyticity in the',& 187 & ' direction (cartesian coordinates)',qphon(1:3)+tol10 188 write(ab_out,'(a)')& 189 & ' -----------------------------------------------------------------------' 190 write(ab_out,*) ch10 191 end if 192 193 else 194 if (iwrite) then 195 write(ab_out,*) ch10 196 write(ab_out,*)' Raman susceptibilities of transverse zone-center phonon modes' 197 write(ab_out,*)' -------------------------------------------------------------' 198 write(ab_out,*) ch10 199 end if 200 end if 201 202 !Examine the degeneracy of each mode. The portability of the echo of the Raman susceptibility 203 !for each degenerate mode is very hard to guarantee. On the contrary, 204 !the scalar reductions of these quantities are OK. 205 ABI_MALLOC(metacharacter,(3*natom)) 206 do imode=1,3*natom 207 ! The degenerate modes are not portable 208 t_degenerate=.false. 209 if(imode>1)then 210 if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true. 211 end if 212 if(imode<3*natom)then 213 if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true. 214 end if 215 metacharacter(imode)=';' 216 if(t_degenerate)metacharacter(imode)='-' 217 end do 218 219 do imode = 1, 3*natom 220 if (iwrite) then 221 write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')' Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)' 222 do idir = 1,3 223 write(ab_out,'(a,4x,3(f16.9,2x))')metacharacter(imode),rsus(imode,idir,:) 224 end do 225 end if 226 ! See R. Caracas and X. Gonze, Thermodynamic Properties of Solids : experiment and modeling, Wiley-VCH, 227 ! Ed. S. Chaplot and R. Mittal and N. Choudhury , chap. 8, pp 291-312. 228 g0=(rsus(imode,1,1)+rsus(imode,2,2)+rsus(imode,3,3))**2*third 229 g1=((rsus(imode,1,2)-rsus(imode,2,1))**2+& 230 & (rsus(imode,1,3)-rsus(imode,3,1))**2+& 231 & (rsus(imode,2,3)-rsus(imode,3,2))**2)*half 232 g2=((rsus(imode,1,2)+rsus(imode,2,1))**2+& 233 & (rsus(imode,1,3)+rsus(imode,3,1))**2+& 234 & (rsus(imode,2,3)+rsus(imode,3,2))**2)*half +& 235 & ((rsus(imode,1,1)-rsus(imode,2,2))**2+& 236 & (rsus(imode,2,2)-rsus(imode,3,3))**2+& 237 & (rsus(imode,3,3)-rsus(imode,1,1))**2)*third 238 if (iwrite) then 239 write(ab_out,'(3(a,f16.9))')' Spherical averages : G0=',g0,' G1=',g1,' G2=',g2 240 write(ab_out,*) 241 end if 242 end do 243 244 ABI_FREE(metacharacter) 245 ABI_FREE(zeff) 246 247 end subroutine ramansus