TABLE OF CONTENTS
ABINIT/ramansus [ Functions ]
NAME
ramansus
FUNCTION
Compute the raman susceptibilities of zone-center phonons
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MVeithen) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 51 subroutine ramansus(d2cart,dchide,dchidt,displ,mpert,& 52 & natom,phfrq,qphon,qphnrm,rsus,ucvol) 53 54 use defs_basis 55 use m_errors 56 use m_profiling_abi 57 58 !This section has been created automatically by the script Abilint (TD). 59 !Do not modify the following lines by hand. 60 #undef ABI_FUNC 61 #define ABI_FUNC 'ramansus' 62 !End of the abilint section 63 64 implicit none 65 66 !Arguments ----------------------------------- 67 !scalars 68 integer,intent(in) :: mpert,natom 69 real(dp),intent(in) :: qphnrm,ucvol 70 !arrays 71 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),dchide(3,3,3) 72 real(dp),intent(in) :: dchidt(natom,3,3,3),displ(2,3*natom,3*natom) 73 real(dp),intent(in) :: phfrq(3*natom) 74 real(dp),intent(inout) :: qphon(3) 75 real(dp),intent(out) :: rsus(3*natom,3,3) 76 77 !Local variables------------------------------- 78 !scalars 79 integer :: analyt,i1,i1dir,i1pert,i2dir,iatom,idir,imode 80 real(dp) :: epsq,fac,g0,g1,g2,qphon2 81 logical :: t_degenerate,iwrite 82 character(len=500) :: message 83 !arrays 84 real(dp) :: dijk_q(3,3) 85 real(dp),allocatable :: zeff(:,:) 86 character(len=1),allocatable :: metacharacter(:) 87 88 ! ********************************************************************* 89 90 iwrite = ab_out > 0 91 92 ABI_ALLOCATE(zeff,(3,natom)) 93 94 rsus(:,:,:) = zero 95 epsq = zero 96 zeff(:,:) = zero 97 dijk_q(:,:) = zero 98 99 !Determine the analyticity of the matrix. 100 analyt=1 101 if(abs(qphnrm)<tol8)analyt=0 102 if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=1 103 104 !In the case the non-analyticity is required : 105 if(analyt == 0) then 106 107 ! Normalize the limiting direction 108 qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2 109 qphon(1)=qphon(1)/sqrt(qphon2) 110 qphon(2)=qphon(2)/sqrt(qphon2) 111 qphon(3)=qphon(3)/sqrt(qphon2) 112 113 ! Get the dielectric constant for the limiting direction 114 epsq= 0._dp 115 do i1dir=1,3 116 do i2dir=1,3 117 epsq=epsq+qphon(i1dir)*qphon(i2dir)*d2cart(1,i1dir,natom+2,i2dir,natom+2) 118 end do 119 end do 120 121 ! Check if epsq > 0 122 if (epsq < tol8) then 123 write(message,'(a,es14.6)')' The value of epsq must be > 0 while it is found to be',epsq 124 MSG_BUG(message) 125 end if 126 127 ! Get the effective charges for the limiting direction 128 do i1dir=1,3 129 do i1pert=1,natom 130 zeff(i1dir,i1pert)=zero 131 do i2dir=1,3 132 zeff(i1dir,i1pert)=zeff(i1dir,i1pert)+qphon(i2dir)*& 133 & d2cart(1,i1dir,i1pert,i2dir,natom+2) 134 end do 135 end do 136 end do 137 138 ! Get the NLO tensor for the limiting direction !$\sum_{k} d_{ijk} \cdot q_k$ 139 140 dijk_q(:,:) = zero 141 do i1dir = 1, 3 142 do i2dir = 1, 3 143 do idir = 1, 3 144 dijk_q(i1dir,i2dir) = dijk_q(i1dir,i2dir) + dchide(i1dir,i2dir,idir)*qphon(idir) 145 end do 146 end do 147 end do 148 149 fac = 16._dp*pi/(ucvol*epsq) 150 do imode = 1, 3*natom 151 do iatom = 1, natom 152 do idir = 1, 3 153 i1=idir + (iatom - 1)*3 154 rsus(imode,:,:) = rsus(imode,:,:) + & 155 & (dchidt(iatom,idir,:,:) - fac*zeff(idir,iatom)*dijk_q(:,:))* & 156 & displ(1,i1,imode) 157 end do ! disp 158 end do ! iatom 159 end do ! imode 160 rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol) 161 162 else 163 do imode = 1, 3*natom 164 do iatom = 1, natom 165 do idir = 1, 3 166 i1=idir + (iatom - 1)*3 167 rsus(imode,:,:) = rsus(imode,:,:) + dchidt(iatom,idir,:,:)*displ(1,i1,imode) 168 end do ! disp 169 end do ! iatom 170 end do ! imode 171 rsus(:,:,:) = rsus(:,:,:)*sqrt(ucvol) 172 end if ! analyt == 0 173 174 if (analyt == 0) then 175 if (iwrite) then 176 write(ab_out,*) ch10 177 write(ab_out, '(a,/,a,3f9.5)' )& 178 & ' Raman susceptibility of zone-center phonons, with non-analyticity in the',& 179 & ' direction (cartesian coordinates)',qphon(1:3)+tol10 180 write(ab_out,'(a)')& 181 & ' -----------------------------------------------------------------------' 182 write(ab_out,*) ch10 183 end if 184 185 else 186 if (iwrite) then 187 write(ab_out,*) ch10 188 write(ab_out,*)' Raman susceptibilities of transverse zone-center phonon modes' 189 write(ab_out,*)' -------------------------------------------------------------' 190 write(ab_out,*) ch10 191 end if 192 end if 193 194 !Examine the degeneracy of each mode. The portability of the echo of the Raman susceptibility 195 !for each degenerate mode is very hard to guarantee. On the contrary, 196 !the scalar reductions of these quantities are OK. 197 ABI_ALLOCATE(metacharacter,(3*natom)) 198 do imode=1,3*natom 199 ! The degenerate modes are not portable 200 t_degenerate=.false. 201 if(imode>1)then 202 if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true. 203 end if 204 if(imode<3*natom)then 205 if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true. 206 end if 207 metacharacter(imode)=';' 208 if(t_degenerate)metacharacter(imode)='-' 209 end do 210 211 do imode = 1, 3*natom 212 if (iwrite) then 213 write(ab_out,'(a4,i3,2x,a2,f7.2,a6)')' Mode',imode,' (',phfrq(imode)*Ha_cmm1,' cm-1)' 214 do idir = 1,3 215 write(ab_out,'(a,4x,3(f16.9,2x))')metacharacter(imode),rsus(imode,idir,:) 216 end do 217 end if 218 ! See R. Caracas and X. Gonze, Thermodynamic Properties of Solids : experiment and modeling, Wiley-VCH, 219 ! Ed. S. Chaplot and R. Mittal and N. Choudhury , chap. 8, pp 291-312. 220 g0=(rsus(imode,1,1)+rsus(imode,2,2)+rsus(imode,3,3))**2*third 221 g1=((rsus(imode,1,2)-rsus(imode,2,1))**2+& 222 & (rsus(imode,1,3)-rsus(imode,3,1))**2+& 223 & (rsus(imode,2,3)-rsus(imode,3,2))**2)*half 224 g2=((rsus(imode,1,2)+rsus(imode,2,1))**2+& 225 & (rsus(imode,1,3)+rsus(imode,3,1))**2+& 226 & (rsus(imode,2,3)+rsus(imode,3,2))**2)*half +& 227 & ((rsus(imode,1,1)-rsus(imode,2,2))**2+& 228 & (rsus(imode,2,2)-rsus(imode,3,3))**2+& 229 & (rsus(imode,3,3)-rsus(imode,1,1))**2)*third 230 if (iwrite) then 231 write(ab_out,'(3(a,f16.9))')' Spherical averages : G0=',g0,' G1=',g1,' G2=',g2 232 write(ab_out,*) 233 end if 234 end do 235 236 ABI_DEALLOCATE(metacharacter) 237 ABI_DEALLOCATE(zeff) 238 239 end subroutine ramansus