TABLE OF CONTENTS
ABINIT/elph2_fanddw [ Functions ]
NAME
elph2_fanddw
FUNCTION
This routine calculates the zero-point motion corrections due to the Fan term or to the DDW term..
COPYRIGHT
Copyright (C) 2011-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors .
INPUTS
dim_eig2nkq=1 if eig2nkq is to be computed displ(2*3*natom*3*natom)=the displacements of atoms in cartesian coordinates. eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=one half second derivatives of the electronic eigenvalues gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) mband= maximum number of bands natom= number of atoms in the unit cell nkpt= number of k-points nsppol= 1 for unpolarized, 2 for spin-polarized option 1 for Fan term, 2 for DDW term phfrq(3*natom)=phonon frequencies (prtvol > 4) if the mode decomposition is to be printed
OUTPUT
eigen_corr(mband*nkpt*nsppol)= T=0 correction to the electronic eigenvalues, due to the Fan term.
PARENTS
respfn
CHILDREN
wrtout
SOURCE
41 #if defined HAVE_CONFIG_H 42 #include "config.h" 43 #endif 44 45 #include "abi_common.h" 46 47 48 subroutine elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_corr,gprimd,mband,natom,nkpt,nsppol,option,phfrq,prtvol) 49 50 use defs_basis 51 use m_errors 52 use m_profiling_abi 53 54 !This section has been created automatically by the script Abilint (TD). 55 !Do not modify the following lines by hand. 56 #undef ABI_FUNC 57 #define ABI_FUNC 'elph2_fanddw' 58 use interfaces_14_hidewrite 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer,intent(in) :: dim_eig2nkq,mband,natom,nkpt,nsppol,option,prtvol 66 67 !arrays 68 real(dp) :: gprimd(3,3) 69 real(dp),intent(in) :: displ(2*3*natom*3*natom) 70 real(dp),intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) 71 real(dp),intent(in) :: phfrq(3*natom) 72 real(dp),intent(out) :: eigen_corr(mband*nkpt*nsppol) 73 74 !Local variables------------------------------- 75 !scalars 76 integer,parameter :: neigs_per_line=6 77 integer :: iatom1,iatom2,idir1,idir2,iband,ikpt,imode,index,isppol, imin, ii 78 real(dp) :: d_at1_dir1_re,d_at1_dir1_im 79 real(dp) :: d_at1_dir2_re,d_at1_dir2_im 80 real(dp) :: d_at2_dir1_re,d_at2_dir1_im 81 real(dp) :: d_at2_dir2_re,d_at2_dir2_im 82 real(dp) :: e2_im,e2_re 83 real(dp), allocatable :: eigen_corr_mode(:) 84 character(len=500) :: message 85 character(len=20) :: eig_format, line_format 86 !arrays 87 real(dp) :: displ2cart(2,3,3),displ2red(2,3,3),tmp_displ2(2,3,3) 88 89 ! ********************************************************************* 90 91 !DEBUG 92 !write(std_out,*)' elph2_fanddw : enter ' 93 !write(std_out,*)' option=',option 94 !ENDDEBUG 95 96 if(option/=1 .and. option/=2)then 97 write(message,'(a,i0)')' The argument option should be 1 or 2, while it is found that option=',option 98 MSG_BUG(message) 99 end if 100 101 !printing options 102 eig_format='f16.8' 103 write(line_format,'(a,i1,a6,a)') '(',neigs_per_line,eig_format,')' 104 105 if (prtvol > 4) then 106 write(message,'(a,a)')ch10,' ================================================================================' 107 call wrtout(ab_out_default,message,'COLL') 108 if (option==1) then 109 write(message,'(a)') ' ---- Begin Fan contributions to eigenvalues renormalization by mode ----' 110 call wrtout(ab_out_default,message,'COLL') 111 else if (option==2) then 112 write(message,'(a)') ' ---- Begin DDW contributions to eigenvalues renormalization by mode ----' 113 call wrtout(ab_out_default,message,'COLL') 114 end if 115 end if 116 117 ABI_ALLOCATE(eigen_corr_mode,(mband*nkpt*nsppol)) 118 119 eigen_corr(:)=zero 120 do imode=1,3*natom 121 eigen_corr_mode(:)=zero 122 123 ! DEBUG 124 ! write(std_out,*)' Contribution of mode ',imode,' with frequency=',phfrq(imode),' and displacements :' 125 ! write(std_out,'(2f14.7)' ) displ(1+2*3*natom*(imode-1):2*3*natom*imode) 126 ! ENDDEBUG 127 128 if (phfrq(imode)>tol6) then 129 do iatom1=1,natom 130 do iatom2=1,natom 131 ! DEBUG 132 ! write(std_out,*)' iatom1,iatom2=',iatom1,iatom2 133 ! ENDDEBUG 134 135 do idir1=1,3 136 do idir2=1,3 137 ! Compute the mean cartesian displacements 138 d_at1_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1)))) 139 d_at1_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1)))) 140 d_at2_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1)))) 141 d_at2_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1)))) 142 143 ! DEBUG 144 ! write(std_out,*)' idir1,idir2=',iatom1,iatom2,idir1,idir2 145 ! write(std_out,'(a,4f12.5)' )' d_at1_dir1 re,d_at2_dir2 re=',d_at1_dir1_re,d_at2_dir2_re 146 ! ENDDEBUG 147 148 if(option==1)then 149 ! Compute the mean displacement correlation at T=0. 150 ! Consistent with Eqs.(7) and (8) of PRB51, 8610 (1995), specialized for the contribution of one q point. 151 ! but generalized to two different atoms. Note that the complex conjugate is taken on the second direction. 152 displ2cart(1,idir1,idir2)=(d_at1_dir1_re*d_at2_dir2_re+ & 153 & d_at1_dir1_im*d_at2_dir2_im )/(two*phfrq(imode)) 154 displ2cart(2,idir1,idir2)=(d_at1_dir1_im*d_at2_dir2_re- & 155 & d_at1_dir1_re*d_at2_dir2_im )/(two*phfrq(imode)) 156 else if(option==2)then 157 ! Compute the mean square displacement correlation of each atom at T=0, and take mean over iatom1 and iatom2. 158 ! See Eqs.(7) and (8) of PRB51, 8610 (1995), specialized for the contribution of one q point. 159 ! Note that the complex conjugate is taken on the second direction. 160 ! Also, note the overall negative sign, to make it opposite to the Fan term. 161 d_at1_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1)))) 162 d_at1_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1)))) 163 d_at2_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1)))) 164 d_at2_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1)))) 165 displ2cart(1,idir1,idir2)=-(d_at1_dir1_re*d_at1_dir2_re+ & 166 & d_at1_dir1_im*d_at1_dir2_im+ & 167 & d_at2_dir1_re*d_at2_dir2_re+ & 168 & d_at2_dir1_im*d_at2_dir2_im )/(four*phfrq(imode)) 169 displ2cart(2,idir1,idir2)=-(d_at1_dir1_im*d_at1_dir2_re- & 170 & d_at1_dir1_re*d_at1_dir2_im+ & 171 & d_at2_dir1_im*d_at2_dir2_re- & 172 & d_at2_dir1_re*d_at2_dir2_im )/(four*phfrq(imode)) 173 end if 174 end do 175 end do 176 ! Switch to reduced coordinates in two steps 177 tmp_displ2(:,:,:)=zero 178 do idir1=1,3 179 do idir2=1,3 180 tmp_displ2(:,:,idir1)=tmp_displ2(:,:,idir1)+displ2cart(:,:,idir2)*gprimd(idir2,idir1) 181 end do 182 end do 183 displ2red(:,:,:)=zero 184 do idir1=1,3 185 do idir2=1,3 186 displ2red(:,idir1,:)=displ2red(:,idir1,:)+tmp_displ2(:,idir2,:)*gprimd(idir2,idir1) 187 end do 188 end do 189 ! Compute the T=0 shift due to this q point 190 do idir1=1,3 191 do idir2=1,3 192 do ikpt=1,nkpt 193 do isppol=1,nsppol 194 do iband=1,mband 195 index=iband+mband*(isppol-1 + nsppol*(ikpt-1)) 196 e2_re=eig2nkq(1,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2) 197 e2_im=eig2nkq(2,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2) 198 eigen_corr(index)=eigen_corr(index)+& 199 & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2) 200 eigen_corr_mode(index)=eigen_corr_mode(index)+& 201 & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2) 202 end do ! band 203 end do ! spin 204 end do ! kpt 205 end do ! dir2 206 end do ! dir1 207 end do ! atom2 208 end do ! atom1 209 end if 210 211 if (prtvol > 4) then 212 ! Print the corrections by mode 213 write(message,'(a,i1)') ' imode= ',imode 214 call wrtout(ab_out_default,message,'COLL') 215 216 do ikpt=1,nkpt 217 do isppol=1,nsppol 218 write(message,'(a,i4,a,i1)')' ikpt= ',ikpt,' ispin= ',isppol 219 call wrtout(ab_out_default,message,'COLL') 220 221 imin = mband * (isppol-1 + nsppol*(ikpt-1)) 222 do ii=0, (mband-1)/neigs_per_line 223 write(message, line_format) (eigen_corr_mode(iband+imin), & 224 & iband = 1 + ii * neigs_per_line, min(mband, (ii+1)*neigs_per_line)) 225 call wrtout(ab_out_default,message,'COLL') 226 end do 227 end do 228 end do 229 end if 230 231 end do ! mode 232 233 if (prtvol > 4) then 234 if (option==1) then 235 write(message,'(a)') ' ---- End Fan contribution to eigenvalues renormalization by mode ----' 236 call wrtout(ab_out_default,message,'COLL') 237 else if (option==2) then 238 write(message,'(a)') ' ---- End DDW contribution to eigenvalues renormalization by mode ----' 239 call wrtout(ab_out_default,message,'COLL') 240 end if 241 write(message,'(a,a)')' ================================================================================', ch10 242 call wrtout(ab_out_default,message,'COLL') 243 end if 244 245 ABI_DEALLOCATE(eigen_corr_mode) 246 247 !DEBUG 248 !write(std_out,*)' elph2_fanddw : exit' 249 !write(std_out,*)' eigen_corr(1)=',eigen_corr(1) 250 !ENDDEBUG 251 252 end subroutine elph2_fanddw