!******************************************************************************* ! * ! Calculate the (optical) dielectric function and joint density of states. * ! Some parts/ideas have been taken from the original codes of Birgit Adolph, * ! however, most of the code has been completely rewritten by J. Furthmueller * ! * !******************************************************************************* PROGRAM OPTICS USE prec USE constant IMPLICIT REAL(q) (A-H,O-Z) INTEGER, PARAMETER :: IU6=6 REAL(q), PARAMETER :: HARTREE=2._q*RYTOEV, THIRD=(1._q/3._q) CHARACTER*72 LINE LOGICAL :: LTET,LJDOS,LDOS,LKRAMERS,LEXTERNAL,LSEARCH INTEGER :: ISYMM COMPLEX(q) :: RESULT,EPS1C,OMEGA INTEGER, ALLOCATABLE :: IDTET(:,:) REAL(q), ALLOCATABLE :: EPSI(:,:,:),EPSR(:,:,:),DOS(:,:), & AMP(:,:,:,:),ETRANS(:,:,:),AMPDOS(:,:,:), & WTKPT(:),PAR(:,:,:,:,:),DOSPAR(:,:,:,:) REAL(q), ALLOCATABLE, TARGET :: FERWE(:,:,:) REAL(q), POINTER :: EIGEN(:,:,:) ! read control flags and customizable parameters CALL READER_OPT(EMIN,EMAX,DE,NEDOS,ISYMM,ISMEAR,SIGMA,LTET, & LJDOS,LDOS,LKRAMERS,LEXTERNAL,LSEARCH,NBCON_USE, & GAMMA,EBOT,ETOP,AZ,FZ,IU6) IF (LEXTERNAL) THEN CALL SCAN_IMEPS(EMAX,DE,NEDOS,IU6) EMIN = 0._q ISPIN = 1 ; NKPTS = 1 ; NTET = 1 ; NBVAL = 1 ; NBCON = 1 ; NBTOT = 2 ELSE ! scan POSCAR, KPOINTS and OPTIC in order to get all dimension parameters CALL SCAN_PARM(IRECL,NKPTS,NTET,ISPIN,NBTOT,NBVAL,NBCON, & LTET,VOL,VOLWGT) CONST = 4._q * PI**2 * HARTREE**3 / VOL ENDIF ! throw away the topmost 15% of conduction bands (usually not well converged)? IF (NBCON_USE==0) NBCON_USE = MAX(INT(0.85_q*NBCON), 1) NBCON_USE = MIN(NBCON_USE, NBCON) ! resulting number of transitions used NTRANS = NBVAL*NBCON_USE IF (ISYMM==1) THEN MAXDIR=1 ELSE IF (ISYMM==2) THEN MAXDIR=2 ELSE IF (ISYMM==3) THEN MAXDIR=3 ELSE MAXDIR=6 ENDIF ! allocate all arrays ALLOCATE( EPSI(NEDOS,ISPIN,MAXDIR), EPSR(NEDOS,ISPIN,MAXDIR), & AMP(NTRANS,NKPTS,ISPIN,6), ETRANS(NTRANS,NKPTS,ISPIN) ) ALLOCATE( IDTET(0:4,NTET), FERWE(NTRANS,NKPTS,ISPIN), WTKPT(NKPTS) ) ALLOCATE( PAR(NTRANS,NKPTS,1,1,ISPIN), DOSPAR(NEDOS,1,1,ISPIN) ) FERWE=0 ; PAR=0 ; DOSPAR=0 ; WTKPT=0 ; IDTET=0 ! during calculation of whatever kind of DOS array FERWE is kept untouched ! and not needed at all, it is just requested as a dummy argument to the ! integration routine with corresponding array size - hence we can reuse it! EIGEN => FERWE(1:NBTOT,:,:) IF (LEXTERNAL) GOTO 100 ! get all input data from files OPTIC and KPOINTS CALL INPUT(IRECL,NTET,NTRANS,NKPTS,NBTOT,NBVAL,NBCON,ISPIN, & LSEARCH,WTKPT(:),AMP,ETRANS,EIGEN(:,:,:),IDTET(:,:)) ! symmetrization IF (ISYMM==1) THEN AMP(:,:,:,1) = (AMP(:,:,:,1) + AMP(:,:,:,2) + AMP(:,:,:,3)) * THIRD ELSE IF (ISYMM==2) THEN AMP(:,:,:,1) = (AMP(:,:,:,1) + AMP(:,:,:,2)) * 0.5_q AMP(:,:,:,2) = AMP(:,:,:,3) ENDIF ! analyze from where transitions come (k points/band pairs versus energy) IF (LSEARCH) THEN CALL SEARCH(NTRANS,NKPTS,NBCON_USE,ETRANS,AMP,ISPIN,MAXDIR,IU0) STOP ENDIF ! DOS and joint DOS: IF (LJDOS.OR.LDOS) THEN write(*,*) 'DOS and joint DOS' call vtime(vpu0,cpu0) ALLOCATE( DOS(NEDOS,ISPIN), AMPDOS(NTRANS,NKPTS,ISPIN) ) AMPDOS=1._q EMINVAL = FLOOR( MINVAL(EIGEN(:,:,:)) - DE ) ! shift for eigenvalues EMAXVAL = MAXVAL(EIGEN(1:NBVAL,:,:)) - EMINVAL ! shifted value of VBM EIGEN = EIGEN - EMINVAL IF (LDOS) THEN IF (LTET) THEN write(*,*) 'Using tetrahedron method' CALL BZINTS(0,FERWE,EIGEN(:,:,:),AMPDOS,WTKPT,NBTOT,NBTOT, & NKPTS,IDTET(:,:),NTET,ISPIN,VOLWGT,EMIN,EMAX, & EPSI(1,1,1),EPSR(1,1,1),NEDOS,EFERMI, & SUMWEI,SUME,IU6,PAR,DOSPAR,NKPTS,1,1,1,0) ELSE if (ismear<0) write(*,*)'Using Lorentzians',SIGMA if (ismear>=0) write(*,*)'Using smearings',SIGMA CALL DENMP(NBTOT,NKPTS,NKPTS,NBTOT,WTKPT,ISPIN, & EMIN,EMAX,ISMEAR,SIGMA,EIGEN(:,:,:),AMPDOS, & NEDOS,EPSI(1,1,1),EPSR(1,1,1)) ENDIF DOS(:,1) = EPSI(:,1,1) IF (ISPIN==2) DOS(:,1) = DOS(:,1) + EPSI(:,ISPIN,1) ELSE DOS=0._q ENDIF IF (LJDOS) THEN IF (LTET) THEN write(*,*) 'Using tetrahedron method' CALL BZINTS(0,FERWE,ETRANS,AMPDOS,WTKPT,NTRANS,NTRANS, & NKPTS,IDTET(:,:),NTET,ISPIN,VOLWGT,EMIN,EMAX, & EPSI(1,1,1),EPSR(1,1,1),NEDOS,EFERMI, & SUMWEI,SUME,IU6,PAR,DOSPAR,NKPTS,1,1,1,0) ELSE if (ismear<0) write(*,*)'Using Lorentzians',SIGMA if (ismear>=0) write(*,*)'Using smearings',SIGMA CALL DENMP(NTRANS,NKPTS,NKPTS,NTRANS,WTKPT,ISPIN, & EMIN,EMAX,ISMEAR,SIGMA,ETRANS,AMPDOS, & NEDOS,EPSI(1,1,1),EPSR(1,1,1)) ENDIF IF (ISPIN==2) EPSI(:,1,1) = EPSI(:,1,1) + EPSI(:,ISPIN,1) ELSE EPSI=0._q ENDIF call vtime(vpu,cpu) vpu=vpu-vpu0 ; cpu=cpu-cpu0 write(*,'(A29,2F10.3)')'Timing DOS/JDOS (cpu/wall): ',vpu,cpu write(*,*) 'output of DOS' OPEN(UNIT = 20, FILE = 'JDOS') WRITE(20,'(A31,F20.16,A13,I8)') '# Joint DOS and DOS shifted by ', & -EMAXVAL,' NEDOS = ',NEDOS DO I=1,NEDOS EACT = EMIN + (I-1)*DE WRITE(20,'(3F24.16)') EACT,EPSI(I,1,1),DOS(I,1) ENDDO CLOSE(20) DEALLOCATE(DOS, AMPDOS) ENDIF 100 EPSI = 0._q ; EPSR = 0._q ; NULLIFY(EIGEN) IF (LEXTERNAL) THEN OPEN(20,FILE='IMEPS',FORM='FORMATTED',STATUS='OLD') GOTO 200 ENDIF 200 CONTINUE dir: DO JDIR = 1, MAXDIR IF (LEXTERNAL) THEN READ(20,'(A72)') LINE DO I=1,NEDOS READ(20,*) EACT,RDUM,EPSI(I,1,JDIR) ENDDO ELSE write(*,*) 'imaginary part of dielectric function' call vtime(vpu0,cpu0) IF (LTET) THEN write(*,*) 'Using tetrahedron method' CALL BZINTS(0,FERWE,ETRANS,AMP(1,1,1,JDIR),WTKPT,NTRANS,NTRANS, & NKPTS,IDTET(:,:),NTET,ISPIN,VOLWGT,EMIN,EMAX, & EPSI(1,1,JDIR),EPSR(1,1,JDIR),NEDOS,EFERMI, & SUMWEI,SUME,IU6,PAR,DOSPAR,NKPTS,1,1,1,0) ELSE if (ismear<0) write(*,*)'Using Lorentzians',SIGMA if (ismear>=0) write(*,*)'Using smearings',SIGMA CALL DENMP(NTRANS,NKPTS,NKPTS,NTRANS,WTKPT,ISPIN, & EMIN,EMAX,ISMEAR,SIGMA,ETRANS,AMP(1,1,1,JDIR), & NEDOS,EPSI(1,1,JDIR),EPSR(1,1,JDIR)) ENDIF call vtime(vpu,cpu) vpu=vpu-vpu0 ; cpu=cpu-cpu0 write(*,'(A29,2F10.3)') 'Timing Im(eps) (cpu/wall): ',vpu,cpu ! rescaling of result to atomic units IF (ISPIN==1) EPSI(:,1,JDIR) = EPSI(:,1,JDIR) * CONST IF (ISPIN==2) EPSI(:,1,JDIR) = (EPSI(:,1,JDIR) + & EPSI(:,ISPIN,JDIR)) * CONST ENDIF ! Real part of epsilon via Kramers-Kronig IF (LKRAMERS) THEN write(*,*) 'real part of dielectric function, be patient ...' call vtime(vpu0,cpu0) DO I=1,NEDOS EACT = EMIN + (I-1)*DE EPS1C= (0._q,0._q) OMEGA= CMPLX( EACT,GAMMA, KIND=q) CALL KRAMERS(OMEGA,EPSI(1,1,JDIR),DE,NEDOS,AZ,FZ,RESULT) EPS1C= RESULT OMEGA= CMPLX(-EACT,GAMMA, KIND=q) CALL KRAMERS(OMEGA,EPSI(1,1,JDIR),DE,NEDOS,AZ,FZ,RESULT) EPS1C= 1._q + EPS1C + RESULT EPSR(I,1,JDIR) = REAL(EPS1C, KIND=q) ENDDO call vtime(vpu,cpu) vpu=vpu-vpu0 ; cpu=cpu-cpu0 write(*,'(A29,2F10.3)') 'Timing Re(eps) (cpu/wall): ',vpu,cpu ELSE EPSR = 0._q ENDIF ENDDO dir IF (LEXTERNAL) CLOSE(20) write(*,*) 'output of dielectric function' OPEN (UNIT = 20, FILE = 'EPS') components: DO JDIR =1, MAXDIR ! sum rule: effective plasma frequency WP=0._q DO I=2,NEDOS EACT = EMIN + (I-1.5_q)*DE WP = WP + EACT*(EPSI(I,1,JDIR) + EPSI(I-1,1,JDIR)) ENDDO WP = 0.5_q*DE*WP * VOL/(TPI*TPI*NBVAL*HARTREE**2) WRITE(20,'(A12,I1,A4,I1,A16,F18.14,A12,I8)') '# Component ',JDIR, & ' of ',MAXDIR,': w_p,eff/w_p = ',SQRT(WP),' NEDOS = ',NEDOS DO I=1,NEDOS EACT = EMIN + (I-1)*DE WRITE(20,'(3F24.16)') EACT,EPSR(I,1,JDIR),EPSI(I,1,JDIR) ENDDO ENDDO components CLOSE(20) END SUBROUTINE SCAN_PARM(IRECL,NKPTS,NTET,ISPIN,NBTOT,NBVAL,NBCON, & LTET,VOL,VOLWGT) USE prec USE constant IMPLICIT NONE ! global variables LOGICAL LTET INTEGER IRECL,NKPTS,NTET,ISPIN,NBTOT,NBVAL,NBCON REAL(q) VOL,VOLWGT ! local variables CHARACTER*80 COMMENT_LINE CHARACTER*1 CSEL INTEGER NKPTS_KPOINTS,NK REAL(q) RIRECL,RISPIN,RNBTOT,RNBVAL,RNBCON,RNKPTS,DUMMY(4) REAL(q) SCALE,A(3,3) write(*,*) 'scanning POSCAR, KPOINTS and OPTIC' ! get the cell volume from file POSCAR (in atomic units!) OPEN(UNIT=20,FILE='POSCAR',FORM='FORMATTED',STATUS='OLD',ERR=100) READ(20,'(A)',END=100,ERR=100) COMMENT_LINE READ(20,*,END=100,ERR=100) SCALE READ(20,*,END=100,ERR=100) A(1,1),A(2,1),A(3,1) READ(20,*,END=100,ERR=100) A(1,2),A(2,2),A(3,2) READ(20,*,END=100,ERR=100) A(1,3),A(2,3),A(3,3) CLOSE(20) A = A*SCALE VOL = A(1,1)*(A(2,2)*A(3,3)-A(3,2)*A(2,3)) + & A(2,1)*(A(3,2)*A(1,3)-A(1,2)*A(3,3)) + & A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) VOL = VOL / (AUTOA*AUTOA*AUTOA) ! now scan file KPOINTS for the number of k-points and tetrahedra OPEN(UNIT=20,FILE='KPOINTS',FORM='FORMATTED',STATUS='OLD',ERR=110) READ(20,'(A)',END=110,ERR=110) COMMENT_LINE READ(20,*,END=110,ERR=110) NKPTS_KPOINTS IF (NKPTS_KPOINTS==0) THEN WRITE(*,*) 'Do not support automatic mode (NKPTS=0). '& & //'Please use explicit list of k-points!' WRITE(*,*) 'To generate such an explicit list launch ' & & //'VASP and copy IBZKPT to KPOINTS.' CLOSE(20) STOP ENDIF READ(20,'(A)',END=110,ERR=110) CSEL DO NK=1,NKPTS_KPOINTS READ(20,*,END=140,ERR=110) DUMMY(1),DUMMY(2),DUMMY(3),DUMMY(4) ENDDO READ(20,'(A)',END=130,ERR=130) CSEL READ(20,*,END=110,ERR=110) NTET,VOLWGT 10 CLOSE(20) ! finally scan file OPTIC OPEN(UNIT=20,FILE='OPTIC',FORM='UNFORMATTED',ACCESS='DIRECT', & STATUS='OLD',RECL=48,ERR=120) READ(20,REC=1,ERR=120) RIRECL, RNBTOT, RNBVAL, RNBCON, RNKPTS, RISPIN CLOSE(20) IRECL =NINT(RIRECL) NBVAL =NINT(RNBVAL) NBCON =NINT(RNBCON) NBTOT =NINT(RNBTOT) NKPTS =NINT(RNKPTS) IF (NKPTS /= NKPTS_KPOINTS) GOTO 140 ISPIN =NINT(RISPIN) GOTO 150 100 WRITE(*,*) 'Error opening or reading file POSCAR.' STOP 110 WRITE(*,*) 'Error opening or reading file KPOINTS.' STOP 120 WRITE(*,*) 'Error opening or reading file OPTIC.' STOP 130 CONTINUE IF (LTET) THEN WRITE(*,*) 'Seems to me there exists no tetrahedra list on file KPOINTS.' WRITE(*,*) 'Please supply one, e.g., launch VASP with automatic k-point' WRITE(*,*) 'generation, set ISMEAR=-4 in file INCAR and use file IBZKPT.' STOP ELSE NTET=0 VOLWGT=0._q ENDIF GOTO 10 140 WRITE(*,*) 'Incompatible number of k-points on file KPOINTS and OPTIC.' STOP 150 CONTINUE RETURN END SUBROUTINE SCAN_IMEPS(EMAX,DE,NEDOS,IU6) USE prec IMPLICIT NONE ! global INTEGER NEDOS,IU6 REAL(q) EMAX,DE ! local CHARACTER*64 LINE INTEGER I REAL(q) DUM1,DUM2 REAL(q), ALLOCATABLE :: E(:) write(*,*) 'scanning IMEPS' OPEN(UNIT=20,FILE='IMEPS',FORM='FORMATTED',STATUS='OLD',ERR=100) READ(20,'(A64,I8)',END=100,ERR=100) LINE,NEDOS ALLOCATE(E(NEDOS)) DO I=1,NEDOS READ(20,*,END=100,ERR=100) E(I),DUM1,DUM2 ENDDO EMAX = E(NEDOS) DE = EMAX / REAL(NEDOS-1, KIND=q) ! consistency check IF (ABS(E(1))>1.E-10_q) GOTO 200 DO I=2,NEDOS IF (ABS(E(I)-E(I-1)-DE)>1.E-10_q) GOTO 200 ENDDO DEALLOCATE(E) CLOSE(20) RETURN 100 WRITE(IU6,*) 'Error opening or reading file IMEPS.' STOP 200 WRITE(IU6,*) 'Internal inconsistency file IMEPS: non-uniform mesh?' STOP END SUBROUTINE INPUT(IRECL,NTET,NTRANS,NKPTS,NBTOT,NBVAL,NBCON,ISPIN, & LSEARCH,WTKPT,AMP,ETRANS,EIGENVAL,IDTET) USE prec USE constant IMPLICIT NONE ! global LOGICAL LSEARCH INTEGER IRECL, NTET, NTRANS, NKPTS, NBTOT, NBVAL, NBCON, ISPIN INTEGER IDTET(0:4,NTET) REAL(q) AMP(NTRANS,NKPTS,ISPIN,6), ETRANS(NTRANS,NKPTS,ISPIN) REAL(q) EIGENVAL(NBTOT,NKPTS,ISPIN), WTKPT(NKPTS) ! local CHARACTER*80 COMMENT_LINE CHARACTER*1 CSEL,CHARAC LOGICAL LQUASI INTEGER I, J, NK, ISP, IDIR, NN, N, NP, IDUM REAL(q) RDUM REAL(q), ALLOCATABLE :: VKPT(:,:),EIGEN(:,:,:),FERWE(:,:,:) REAL(q), ALLOCATABLE :: VKPT_KPOINTS(:,:),WTKPT_KPOINTS(:) REAL(q), ALLOCATABLE :: QUASI(:,:,:) COMPLEX(q), ALLOCATABLE :: NABIJ(:,:,:,:,:) ALLOCATE( VKPT(3,NKPTS), & EIGEN(NBTOT,NKPTS,ISPIN), FERWE(NBTOT,NKPTS,ISPIN), & QUASI(NBTOT,NKPTS,ISPIN) ) ALLOCATE( VKPT_KPOINTS(3,NKPTS), WTKPT_KPOINTS(NKPTS) ) ALLOCATE( NABIJ(NBVAL,3,ISPIN,NKPTS,NBCON) ) write(*,*) 'reading KPOINTS and OPTIC' ! tetrahedron list OPEN(UNIT=20, FILE='KPOINTS', FORM='FORMATTED', STATUS='OLD') READ(20,'(A)') COMMENT_LINE READ(20,*) IDUM READ(20,'(A)') CSEL DO NK=1, NKPTS READ(20,*) (VKPT_KPOINTS(I,NK),I=1,3), WTKPT_KPOINTS(NK) ENDDO READ(20,'(A)',ERR=10,END=10) CHARAC READ(20,*) IDUM,RDUM DO I=1, NTET READ(20, *) (IDTET(J,I), J=0,4) ENDDO 10 CLOSE(20) ! transition matrix elements and LDA eigenvalues OPEN(UNIT=20,FILE='OPTIC',FORM='UNFORMATTED', & ACCESS='DIRECT',RECL=IRECL,STATUS='OLD') DO NK = 1, NKPTS READ(20,REC=NK+1) (VKPT(I,NK),I=1,3), WTKPT(NK) DO ISP = 1, ISPIN READ(20,REC=ISP*NKPTS+NK+1) & (EIGEN(N,NK,ISP),FERWE(N,NK,ISP),N=1,NBVAL) DO NN = 1, NBCON N = NBTOT - NBCON + NN DO IDIR = 1, 3 READ(20,REC=(1+ISPIN)*NKPTS+1+IDIR+3*(ISP-1) & +3*ISPIN*(NK-1)+3*ISPIN*NKPTS*(NN-1)) & EIGEN(N,NK,ISP), FERWE(N,NK,ISP), & (NABIJ(NP,IDIR,ISP,NK,NN),NP=1,NBVAL) ENDDO ENDDO ENDDO ENDDO CLOSE(20) ! make (apparently or really) degenerate eigenvalues to be really degenerate: DO ISP = 1, ISPIN DO NK = 1, NKPTS DO N = 2, NBTOT IF (ABS(EIGEN(N,NK,ISP)-EIGEN(N-1,NK,ISP))<1.E-4_q) & EIGEN(N,NK,ISP)=EIGEN(N-1,NK,ISP) ENDDO ENDDO ENDDO ! quasi-particle shifts from GW calculations if available: QUASI = 0._q INQUIRE(FILE='QPSHIFT',EXIST=LQUASI) IF (LQUASI) THEN OPEN(UNIT=20, FILE='QPSHIFT', FORM='FORMATTED', STATUS='OLD') READ(20,'(A)',ERR=100,END=100) COMMENT_LINE READ(20,'(A)',ERR=100,END=100) COMMENT_LINE READ(20,'(A)',ERR=100,END=100) COMMENT_LINE READ(20,'(A)',ERR=100,END=100) COMMENT_LINE READ(20,'(A)',ERR=100,END=100) COMMENT_LINE READ(20,*,ERR=100,END=100) IDUM,NK,NN IF (NK/=NKPTS) GOTO 110 IF (NN/=NBTOT) GOTO 120 DO NK=1,NKPTS READ(20,*,ERR=100,END=100) READ(20,'(A)',ERR=100,END=100) COMMENT_LINE DO N=1,NBTOT READ(20,*,ERR=100,END=100) IDUM,(QUASI(N,NK,ISP),ISP=1,ISPIN) ENDDO ENDDO CLOSE(20) GOTO 200 100 CONTINUE WRITE(*,*) 'Error opening or reading file QPSHIFT.' CLOSE(20) STOP 110 CONTINUE WRITE(*,*) 'Incompatible number of k points on file QPSHIFT.' CLOSE(20) STOP 120 CONTINUE WRITE(*,*) 'Incompatible number of bands on file QPSHIFT.' CLOSE(20) STOP 200 CONTINUE ENDIF ! quasi-particle eigenvalues entering the delta function in DOS and JDOS EIGENVAL = EIGEN + QUASI CALL POST_INP(NKPTS,NBVAL,NBCON,NBTOT,NTRANS,ISPIN, & LSEARCH,EIGEN,FERWE,NABIJ,AMP,ETRANS,EIGENVAL) DEALLOCATE(VKPT,EIGEN,FERWE,QUASI,VKPT_KPOINTS,WTKPT_KPOINTS,NABIJ) RETURN END SUBROUTINE POST_INP(NKPTS,NBVAL,NBCON,NBTOT,NTRANS,ISPIN, & LSEARCH,EIGEN,FERWE,NABIJ,AMP,ETRANS,EIGENVAL) USE prec USE constant IMPLICIT NONE ! global LOGICAL LSEARCH INTEGER NTRANS, NBTOT, NBVAL, NBCON, NKPTS, ISPIN REAL(q) EIGEN(NBTOT,NKPTS,ISPIN), FERWE(NBTOT,NKPTS,ISPIN), & ETRANS(NTRANS,NKPTS,ISPIN), EIGENVAL(NBTOT,NKPTS,ISPIN), & AMP(NTRANS,NKPTS,ISPIN,6) COMPLEX(q) NABIJ(NBVAL,3,ISPIN,NKPTS,NBCON) ! local INTEGER I, J, IK, ISP, IDIR1, IDIR2, JDIR, NC, NV, NBCON_USE write(*,*) 'set up transition energies and amplitudes' AMP = 0._q ETRANS= 0._q ! maybe throw away the topmost conduction bands (usually not well converged) NBCON_USE = NTRANS / NBVAL DO ISP = 1, ISPIN DO IK = 1, NKPTS DO NV = 1, NBVAL DO NC = 1, NBCON_USE J = (NV-1)*NBCON_USE+ NC DO IDIR2 = 1, 3 DO IDIR1 = IDIR2, 3 ! the following looks very cryptic but it simply performs an index mapping ! 11 -> 1 , 22 -> 2 , 33 -> 3, 21/12 -> 4, 31/13 ->6 and 32/23 -> 6 JDIR=IDIR2+((IDIR1-IDIR2+1)/2)*(IDIR1+1) AMP(J,IK,ISP,JDIR) = AMP(J,IK,ISP,JDIR) + & (NABIJ(NV,IDIR2,ISP,IK,NC)*CONJG(NABIJ(NV,IDIR1,ISP,IK,NC))) ENDDO ENDDO IF (LSEARCH) THEN ! in the search mode we just want to output in kinetic-energy ! equivalents (eV units - requires a scaling with RYTOEV because from VASP ! we get

in 1/a.u., so we need 1/2**HARTREE = *RYTOEV) DO JDIR = 1, 6 AMP(J,IK,ISP,JDIR) = AMP(J,IK,ISP,JDIR) * RYTOEV ENDDO ELSE ! scaling of the matrix elements with the LDA energy denominator ETRANS(J,IK,ISP) = EIGEN(NC+NBTOT-NBCON,IK,ISP) - EIGEN(NV,IK,ISP) DO JDIR = 1, 6 AMP(J,IK,ISP,JDIR) = AMP(J,IK,ISP,JDIR) / & (ETRANS(J,IK,ISP)*ETRANS(J,IK,ISP)) ENDDO ENDIF ! set up the correct transition energies entering the delta function ETRANS(J,IK,ISP) = EIGENVAL(NC+NBTOT-NBCON,IK,ISP) - EIGENVAL(NV,IK,ISP) ENDDO ENDDO ENDDO ENDDO RETURN END SUBROUTINE KRAMERS(X,AM,DE,N,A,FZ,FC) USE prec USE constant ! Computes the integral of dY F(Y)/(Y-X) from A to N*DE ! Special highly optimized version of Birgits code for uniform meshes ! The mesh spacing is given by DE ! The result is FC ! The integrand is in AM(:) ! its value at A is FZ ! is must hold AIMAG(X)>0 IMPLICIT NONE REAL(q), PARAMETER :: CO=1._q, CH=-0.5_q, CT=1._q/3._q, CQ=-0.25_q REAL(q), PARAMETER :: ZCRIT1=2.E-3_q, ZCRIT2=1.E-2_q ! global INTEGER :: N REAL(q) :: AM(N),A,FZ,DE COMPLEX(q) :: X,FC ! local INTEGER , SAVE :: J,N1MIN,N2MIN,N1MAX,N2MAX REAL(q) , SAVE :: ZABS,T COMPLEX(q), SAVE :: FC2,B,C,D,F,H,Z IF (ABS(A)<1.E-10_q) GOTO 111 F=A-X B=LOG(1._q-A/F) FC=FZ*(-1._q+LOG((-X)/F)*(X/A)) ! code for J=1 (taken out of J-loop below since first point may be non-uniform) H=-X D=DE-X Z=D/H C=LOG(Z) FC=FC+AM(1)*(C*D/DE+B*F/A) B=C F=H FC2=(0._q,0._q) ! Since the evaluation of LOG(X) is rather time consuming it is more efficient ! to replace the LOG function call by a (4th/3rd order) Taylor series where the ! argument of the LOG function is close to one. Therefore, the remaining loop ! J=2,N-1 is divided into appropriate regions where the LOG function has to be ! evaluated exactly or where it can be replaced by a simple 4th or 3rd order ! Taylor series. The regions are defined by the parameters ZCRIT1 and ZCRIT2. ! With the current setting one should expect an accuray of the order 1E-8 or ! better. The speedup can be significant (e.g. a factor 8 on an IBM RS6000). T =REAL(X, KIND=q)/DE + 1._q N1MIN= FLOOR(T-1._q/ZCRIT1) N2MIN= FLOOR(T-1._q/ZCRIT2) N2MAX=CEILING(T+1._q/ZCRIT2) N1MAX=CEILING(T+1._q/ZCRIT1) ! region 1: large negative values H, D/H very close to 1, 3rd order series DO J=2,MIN(N1MIN,N-1) H=D D=D+DE Z=DE/H C=Z*(CO+Z*(CH+Z*CT)) FC2=FC2+AM(J)*(C*D-B*F) B=C F=H ENDDO ! region 2: intermediate negative values H, D/H close to 1, 4th order series DO J=MAX(N1MIN+1,2),MIN(N2MIN,N-1) H=D D=D+DE Z=DE/H C=Z*(CO+Z*(CH+Z*(CT+Z*CQ))) FC2=FC2+AM(J)*(C*D-B*F) B=C F=H ENDDO ! region 3: small values of H, D/H rather distant to 1, use exact LOG function DO J=MAX(N2MIN+1,2),MIN(N2MAX-1,N-1) H=D D=D+DE Z=D/H C=LOG(Z) FC2=FC2+AM(J)*(C*D-B*F) B=C F=H ENDDO ! region 4: intermediate positive values H, D/H close to 1, 4th order series DO J=MAX(N2MAX,2),MIN(N1MAX-1,N-1) H=D D=D+DE Z=DE/H C=Z*(CO+Z*(CH+Z*(CT+Z*CQ))) FC2=FC2+AM(J)*(C*D-B*F) B=C F=H ENDDO ! region 5: large positive values H, D/H very close to 1, 3rd order series DO J=MAX(N1MAX,2),N-1 H=D D=D+DE Z=DE/H C=Z*(CO+Z*(CH+Z*CT)) FC2=FC2+AM(J)*(C*D-B*F) B=C F=H ENDDO ! add last point, perform global sum and divide by PI FC=(FC+FC2/DE+AM(N)*(1._q-B*F/DE))/PI RETURN 111 WRITE(*,112) 112 FORMAT(' Error subroutine KRAMERS: A too small.') STOP END SUBROUTINE READER_OPT(EMIN,EMAX,DE,NEDOS,ISYMM,ISMEAR,SIGMA, & LTET,LJDOS,LDOS,LKRAMERS,LEXTERNAL,LSEARCH, & NBCON_USE,GAMMA,EBOT,ETOP,AZ,FZ,IU0) USE prec IMPLICIT NONE ! global LOGICAL LTET,LJDOS,LDOS,LKRAMERS,LEXTERNAL,LSEARCH INTEGER NEDOS,ISYMM,ISMEAR,NBCON_USE,IU0 REAL(q) SIGMA,EMIN,EMAX,DE,GAMMA,EBOT,ETOP,AZ,FZ ! local CHARACTER*1 CHARAC LOGICAL LOPEN,LDUM INTEGER N,IERR,IU5,NXTFRU,IDUM REAL(q) RDUM COMPLEX(q) CDUM LOPEN = .TRUE. IU5 = NXTFRU() ! Minimum energy EMIN = 0._q CALL RDATAB(LOPEN,'OPTCTR',IU5,'OMMIN','=','#',';','F', & & IDUM,EMIN,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''OMMIN'' from file OPTCTR.' GOTO 150 ENDIF ! Maximum energy EMAX = 20._q CALL RDATAB(LOPEN,'OPTCTR',IU5,'OMMAX','=','#',';','F', & & IDUM,EMAX,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''OMMAX'' from file OPTCTR.' GOTO 150 ENDIF ! number of energy intervals NEDOS = 4000 CALL RDATAB(LOPEN,'OPTCTR',IU5,'NEDOS','=','#',';','I', & & NEDOS,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''NEDOS'' from file OPTCTR.' GOTO 150 ENDIF NEDOS = NEDOS+1 DE = (EMAX-EMIN)/REAL(NEDOS-1,KIND=q) ! maximum number of conduction bands, NBCON_USE = 0 means INT(0.85*NBCON) NBCON_USE = 0 CALL RDATAB(LOPEN,'OPTCTR',IU5,'NBCON','=','#',';','I', & & NBCON_USE,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''NBCON'' from file OPTCTR.' GOTO 150 ENDIF ! re-symmetrization of the matrix elements (default: 'cubic') ISYMM=1 CALL RDATAB(LOPEN,'OPTCTR',IU5,'ISYMM','=','#',';','I', & & ISYMM,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''ISYMM'' from file OPTCTR.' GOTO 150 ENDIF ! smearing parameter ISMEAR (default is Lorentzian smearing!) ISMEAR = -1 CALL RDATAB(LOPEN,'OPTCTR',IU5,'ISMEAR','=','#',';','I', & & ISMEAR,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''ISMEAR'' from file OPTCTR.' GOTO 150 ENDIF IF ((ISMEAR.EQ.-2).OR.(ISMEAR.EQ.-3).OR.(ISMEAR<-5)) ISMEAR=-1 IF (ISMEAR.EQ.-5) ISMEAR=-4 ! Smearing parameter SIGMA (default is 0.1 eV) SIGMA = 0.1_q CALL RDATAB(LOPEN,'OPTCTR',IU5,'SIGMA','=','#',';','F', & & IDUM,SIGMA,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''SIGMA'' from file OPTCTR.' GOTO 150 ENDIF ! logical control flag LTET, LJDOS, LDOS, LKRAMERS, LEXTERNAL, and LSEARCH LTET = .TRUE. CALL RDATAB(LOPEN,'OPTCTR',IU5,'LTET','=','#',';','L', & & IDUM,RDUM,CDUM,LTET,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LTET'' from file OPTCTR.' GOTO 150 ENDIF ! ISMEAR = -4 (or -5 matched back to -4 above) is a synonym for LTET=.TRUE. ! and takes precedence over any setting of LTET by logical flag LTET ! IF (ISMEAR.EQ.-4) THEN IF (.NOT.LTET) & WRITE(IU0,*)'Warning! ISMEAR = -4/-5 enforces LTET=.TRUE.' IF (.NOT.LTET) & WRITE(IU0,*)' LTET changed from .FALSE. to .TRUE.' IF (.NOT.LTET) & WRITE(IU0,*)' I will use the tetrahedron method now!' LTET=.TRUE. ENDIF LJDOS = .TRUE. CALL RDATAB(LOPEN,'OPTCTR',IU5,'LJDOS','=','#',';','L', & & IDUM,RDUM,CDUM,LJDOS,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LJDOS'' from file OPTCTR.' GOTO 150 ENDIF LDOS = LJDOS CALL RDATAB(LOPEN,'OPTCTR',IU5,'LDOS','=','#',';','L', & & IDUM,RDUM,CDUM,LDOS,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LDOS'' from file OPTCTR.' GOTO 150 ENDIF LKRAMERS = .TRUE. CALL RDATAB(LOPEN,'OPTCTR',IU5,'LKRAMERS','=','#',';','L', & & IDUM,RDUM,CDUM,LKRAMERS,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LKRAMERS'' from file OPTCTR.' GOTO 150 ENDIF LEXTERNAL = .FALSE. CALL RDATAB(LOPEN,'OPTCTR',IU5,'LEXTERNAL','=','#',';','L', & & IDUM,RDUM,CDUM,LEXTERNAL,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LEXTERNAL'' from file OPTCTR.' GOTO 150 ENDIF IF (LEXTERNAL) LKRAMERS = .TRUE. LSEARCH = .FALSE. CALL RDATAB(LOPEN,'OPTCTR',IU5,'LSEARCH','=','#',';','L', & & IDUM,RDUM,CDUM,LSEARCH,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''LSEARCH'' from file OPTCTR.' GOTO 150 ENDIF ! imaginary shift for Kramers-Kronig GAMMA=.002_q CALL RDATAB(LOPEN,'OPTCTR',IU5,'GAMMA','=','#',';','F', & & IDUM,GAMMA,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''GAMMA'' from file OPTCTR.' GOTO 150 ENDIF IF (GAMMA<1.E-8_q) THEN WRITE(IU0,*)'Warning reading item ''GAMMA'' from file OPTCTR.' WRITE(IU0,*)'Supplied value too small, ''GAMMA'' reset to 1.E-8.' GAMMA=1.E-8_q ENDIF ! lower integration bound on real axis for Kramers-Kronig ETOP = 0.05_q EBOT = 2._q CALL RDATAB(LOPEN,'OPTCTR',IU5,'ETOP','=','#',';','F', & & IDUM,ETOP,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''ETOP'' from file OPTCTR.' GOTO 150 ENDIF CALL RDATAB(LOPEN,'OPTCTR',IU5,'EBOT','=','#',';','F', & & IDUM,EBOT,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''EBOT'' from file OPTCTR.' GOTO 150 ENDIF AZ = 0.5_q*(EBOT-ETOP)/REAL(NEDOS-1, KIND=q) CALL RDATAB(LOPEN,'OPTCTR',IU5,'AZ','=','#',';','F', & & IDUM,AZ,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''AZ'' from file OPTCTR.' GOTO 150 ENDIF ! function value Im(eps) at AZ (assume semiconductor, i.e. a gap => always 0 !) FZ = 0._q RETURN 150 CONTINUE STOP END SUBROUTINE SEARCH(NTRANS,NKPTS,NBCON_USE,ETRANS,AMP,ISPIN,MAXDIR,IU0) USE prec IMPLICIT NONE ! global INTEGER NTRANS, NKPTS, NBCON_USE, ISPIN, MAXDIR, IU0 REAL(q) ETRANS(NTRANS,NKPTS,ISPIN), AMP(NTRANS,NKPTS,ISPIN,6) ! local CHARACTER*1 CHARAC CHARACTER*80 COMMENT_LINE LOGICAL LDUM, LOPEN INTEGER N, IK, ISP, IDIR, NC, NV, IDUM, IERR, IU5, NXTFRU, KMIN, KMAX REAL(q) EMINSEARCH, EMAXSEARCH, AMPMIN, RDUM COMPLEX(q) CDUM REAL(q), ALLOCATABLE :: VKPT(:,:),WTKPT(:) LOPEN = .TRUE. IU5 = NXTFRU() ! Minimum energy for search EMINSEARCH = -1.E30_q CALL RDATAB(LOPEN,'OPTCTR',IU5,'EMINSEARCH','=','#',';','F', & & IDUM,EMINSEARCH,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''EMINSEARCH'' from file OPTCTR.' GOTO 150 ENDIF ! Maximum energy for search EMAXSEARCH = 1.E30_q CALL RDATAB(LOPEN,'OPTCTR',IU5,'EMAXSEARCH','=','#',';','F', & & IDUM,EMAXSEARCH,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''EMAXSEARCH'' from file OPTCTR.' GOTO 150 ENDIF ! Minimum amplitude for printout AMPMIN = 0._q CALL RDATAB(LOPEN,'OPTCTR',IU5,'AMPMIN','=','#',';','F', & & IDUM,AMPMIN,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''AMPMIN'' from file OPTCTR.' GOTO 150 ENDIF ! First k point KMIN=1 CALL RDATAB(LOPEN,'OPTCTR',IU5,'KMIN','=','#',';','I', & & KMIN,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''KMIN'' from file OPTCTR.' GOTO 150 ENDIF ! Last k point KMAX=NKPTS CALL RDATAB(LOPEN,'OPTCTR',IU5,'KMAX','=','#',';','I', & & KMAX,RDUM,CDUM,LDUM,CHARAC,N,1,IERR) IF (((IERR/=0).AND.(IERR/=3)).OR. & & ((IERR==0).AND.(N<1))) THEN IF (IU0>=0) & WRITE(IU0,*)'Error reading item ''KMAX'' from file OPTCTR.' GOTO 150 ENDIF ALLOCATE( VKPT(3,NKPTS), WTKPT(NKPTS) ) ! k point list OPEN(UNIT=20, FILE='KPOINTS', FORM='FORMATTED', STATUS='OLD') READ(20,'(A)') COMMENT_LINE READ(20,*) IDUM READ(20,'(A)') CHARAC DO IK=1, NKPTS READ(20,*) (VKPT(N,IK),N=1,3), WTKPT(IK) ENDDO DO ISP=1,ISPIN DO IK=KMIN,KMAX DO N=1,NTRANS NV=((N-1)/NBCON_USE)+1 NC=MOD(N-1,NBCON_USE)+1 DO IDIR=1,MAXDIR IF ((ETRANS(N,IK,ISP)>=EMINSEARCH).AND. & (ETRANS(N,IK,ISP)<=EMAXSEARCH).AND. & (AMP(N,IK,ISP,IDIR)>=AMPMIN)) THEN WRITE(*,*) WRITE(*,'(2(A,I6),2(A,I2),A,I8,A)') & 'NV =',NV,' NC =',NC,' ISP =',ISP,' IDIR =',IDIR,' IK =',IK,':' WRITE(*,'(2(A,F8.5),A,F7.4,A,F7.4,A,F7.4)') & 'E =',ETRANS(N,IK,ISP),' AMP =',AMP(N,IK,ISP,IDIR), & ' K =',VKPT(1,IK),',',VKPT(2,IK),',',VKPT(3,IK) ENDIF ENDDO ENDDO ENDDO ENDDO DEALLOCATE(VKPT,WTKPT) RETURN 150 CONTINUE STOP END