C (C) Copr. 10/1999 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wave_scattering@yahoo.com C C Last update: October 1999 C C This is LDOSL1D, a computer code to calculate the local C density of states for a 1D toy model. This code is based C on the article C A. Moroz, Minima and maxima of the local C density of states for one-dimensional periodic systems, C Europhys. Lett. 46, 419-424 (1999). C C This code is made available on the following conditions: C C - You can use and modify it, as long as acknowledgment is given C to the author (that's me). I'd appreciate a reprint or e-copy C of any paper where my code has been used. C C - You can not sell it or use it in any way to get money out of it C (save for the research grants you might get thanks to it) C C - I'm making this code available to the best of my knowledge. C It has been tested as throughly as possible. C C - Disclaimer. The usual stuff. I claim no responsibility for any C misuse or damage arising from its use. If this program erases your C hard disk, makes your computer cry out for Linux, has you fired, C makes you lose your boy/girlfriend, melts your TV set, spends out C your lifesaving or forces you to get to decaf coffee, it's your C problem. For your tranquility, no such side effects has been C reported ... yet. C C Improvements, comments, additions: all welcome at C C wave_scattering@yahoo.com C PROGRAM ldosl1d C---------------------------------------------------------------------- C PROGRAM TO CALCULATE THE LOCAL DOS OF PHOTONS/ELECTRONS C ON A 1D LATTICE NEAR A BAND-EDGE C C f77 -g -check_bounds ldosl1d.f -o rnldosl1d C C Produces files: C ldsedj.dat where j=1,2, ... 5 - contain the LDOS asymptotic near a C band edge for 5 equidistant points C within a unit cell. C File ldsed1.dat is for the point A, the file ldsed5.dat for the point E C C In the 1st column ... a distance from the band edge under considerations C the 2nd column ... the LDOS C the 3rd column ... the ratio of the log of the LDOS to the distance C from the band edge C C Note: C The program parameters correspond to a 1D structure discussed in the article. C The structure consists of two alternating layers with their respective width C ratio $3/2$ and their respective dielectric constants $1$ and $12$. C Unit cells are centered around the layer with the higher dielectric constant C (called scatterer in the program). For this configuration, the edges of the C first gap are at frequencies C 1.047929335868379 and 1.9879393737078342, respectively. C The edges of the second gap are then at frequencies C 2.6626891316403768 and 3.750103714882322, respectively. C You can run the program for other choice C of parameters, however, then you supply yourself the position of band edges. C C Parameters which you can modify: C YNUPP ... Is the current band edge the upper band edge? C EDW ... the position of a band edge under considerations C ENW ... the first 7 digits of the band edge C XENW ... the last digits of the band edge starting from the 7th decimal place C XSHIFT ... using this parameter you can offset positions of the C 5 equidistant points A-E within the unit cell. C Restriction: 0.GE.XSHIFT.LT. ALA/2.d0, where ALA is the lattice C boundary and constant (set to unity). C If XSHIFT=0, then the point A is at a unit cell boundary and C the point E is exactly at the unit cell center C (in the middle of the scattererr). C EPS0 ... ``host" permittivity C EPS1 ... ``scatterer" permittivity C RMUF ... the scatterer width C C Parameters which is not recommended to alter: C ALA ... the lattice constant (set to unity) C NSTEP ... number of steps on frequency interval C NX ... number of steps on direct lattice C C Arrays: C DOS ... contains the LDOS at a given position within the unit cell C and for a given frequency C---------------------------------------------------------------------- implicit none real*8 ENW,XENW,EDW,XSHIFT real*8 ALA,ADELO,RMUF,EPS0,EPS1 integer*4 NOUT,NSTEP,NX character*1 ynupp C ::: number of the output unit PARAMETER (NOUT=30) C ::: the band edge under considerations PARAMETER (EDW=2.6626891316403768d0) C ::: the first 7 digits of the band edge PARAMETER (ENW=2.662689d0) C ::: the last digits of the band edge starting from the 7th decimal place PARAMETER (XENW=1.316403768d0) C ::: Is the current band edge the upper BAND (not gap) edge? PARAMETER (YNUPP='y') C ::: lattice constant: PARAMETER (ala=1.d0) C ::: muffin-tin radius (rmuf < ala/2.) PARAMETER (rmuf=0.2d0) C ::: BACKGROUND PERMITTIVITY PARAMETER (EPS0=1.D0) C ::: SCATTERER PERMITTIVITY PARAMETER (EPS1=12.D0) C ::: number of steps on direct lattice -choose an even number: PARAMETER (NX=8) C ::: a shift of a step on direct lattice (xshift < ala/2.): PARAMETER (XSHIFT=0.000d0*ala/nx) * real*8 taui(2),rl(2),cpsh(2),spsh(2),tpsh(2) real*8 ff(2),dos(2,0:nx),ddos(0:nx),xcl(2) real*8 pi,p,omega(2),sg0,sg1,phi real*8 xx,xd,xn,xna,xnb,xf,ttf real*8 rl1,rdl1,rl2,rdl2,gamma1,gamma2 real*8 x,xx0,xx1,delx,xnw,rmg integer igf,jx,id,l,idm integer iedge DATA PI/3.141592653589793d0/ if (rmuf .ge. ala/2) then write(6,*) 'Warning: rmuf.ge.ala/2 !!!' stop end if if (xshift .ge. ala/2) then write(6,*) 'Warning: xshift.ge.ala/2 !!!' stop end if WRITE(6,*)'XSHIFT=',XSHIFT OPEN(UNIT=NOUT+1,FILE='ldsed1.dat') rewind(NOUT+1) OPEN(UNIT=NOUT+2,FILE='ldsed2.dat') rewind(NOUT+2) OPEN(UNIT=NOUT+3,FILE='ldsed3.dat') rewind(NOUT+3) OPEN(UNIT=NOUT+4,FILE='ldsed4.dat') rewind(NOUT+4) OPEN(UNIT=NOUT+5,FILE='ldsed5.dat') rewind(NOUT+5) * p=sqrt(eps1/eps0) * *oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo if((ynupp.eq.'y').or.(ynupp.eq.'Y')) then iedge=0 else iedge=1 end if * Scanning over frequency interval: * The exponent below cannot be larger than 7. Otherwise ifix (int) * assignement gives wrong result since integers are only * written on up to 9 digits do 200 id=0,7 idm=1+mod(id,2) xnw=ifix(XENW*10**id)+iedge c write(6,*)'XNW=',xnw omega(idm)=enw + xnw/dble(10.d0**(id+7)) c write(6,*)'XNW=',10.d0**(id+7) c write(6,*)'ID=',id write(nout,*)'omega=', omega(idm) sg0=omega(idm)*sqrt(eps0)*rmuf sg1=omega(idm)*sqrt(eps1)*rmuf * * * xcl(1)= sqrt( cos(sg1)**2 + p**2 *sin(sg1)**2) xcl(2)= sqrt( sin(sg1)**2 + p**2 *cos(sg1)**2) cpsh(1)= ( cos(sg0)*cos(sg1) + p*sin(sg0)*sin(sg1))/xcl(1) spsh(1)= (-sin(sg0)*cos(sg1) + p*cos(sg0)*sin(sg1))/xcl(1) * cpsh(2)= (sin(sg0)*sin(sg1) + p*cos(sg0)*cos(sg1))/xcl(2) spsh(2)= (cos(sg0)*sin(sg1) - p*sin(sg0)*cos(sg1))/xcl(2) * tpsh(1)= spsh(1)/cpsh(1) tpsh(2)= spsh(2)/cpsh(2) * phi=omega(idm)*sqrt(eps0)*ala ff(1) = cos(phi) - tpsh(1)*sin(phi) ff(2) = cos(phi) - tpsh(2)*sin(phi) * xd = cos(psh(1) - psh(2)) * xna= - cos(psh(1) + psh(2)) * xnb= sin(psh(1) + psh(2)) xd = cpsh(1)*cpsh(2)+ spsh(1)*spsh(2) xna= -cpsh(1)*cpsh(2)+ spsh(1)*spsh(2) xnb= cpsh(1)*spsh(2)+ spsh(1)*cpsh(2) xn = -cos(phi)*xna-sin(phi)*xnb xf=xn/xd * logarithmic derivatives at zone boundary: xx = omega(idm)*sqrt(eps0)*ala/2.d0 rl1 = cos(xx)*cpsh(1) - sin(xx)*spsh(1) rdl1 = - sin(xx)*cpsh(1) - cos(xx)*spsh(1) rl2 = sin(xx)*cpsh(2) + cos(xx)*spsh(2) rdl2 = cos(xx)*cpsh(2) - sin(xx)*spsh(2) gamma1=omega(idm)*sqrt(eps0)*rdl1/rl1 gamma2=omega(idm)*sqrt(eps0)*rdl2/rl2 if ( (gamma1-gamma2).le.0. ) then igf=1 else igf=-1 end if if (abs(xf).lt. 1.) then ttf = (1.d0+tpsh(1)*tpsh(2))*sqrt(1.d0-xf**2) taui(1)= omega(idm)*sqrt(eps0)*tpsh(1)*(xf-ff(2))/ttf taui(2)= omega(idm)*sqrt(eps0)*tpsh(2)*(xf-ff(1))/ttf else taui(1)= 0.d0 taui(2)= 0.d0 end if *---------------------------------------------------------------------- * Scanning over elementary cell: delx=ala/dble(nx) * using that the LDOS is symmetric function within the unit cell: * XSHIFT below enables to shift within DELX do 100 jx=0,nx/2 x= -ala/2.d0 + dble(jx)*delx + xshift xx0 = omega(idm)*sqrt(eps0)*abs(x) xx1 = omega(idm)*sqrt(eps1)*abs(x) * Radial solutions if (abs(x).le.rmuf) then * within the muffin-tin radius: rl(1) = p*cos(xx1)/xcl(1) rl(2) = p*sin(xx1)/xcl(2) else * outside muffin-tin radius: rl(1) = (cos(xx0)*cpsh(1) - sin(xx0)*spsh(1)) rl(2) = (sin(xx0)*cpsh(2) + cos(xx0)*spsh(2)) end if * Calculation of the local DOS dos(idm,jx) = 0.d0 do 10 l=1,2 dos(idm,jx) = dos(idm,jx)+ taui(l)* & (rl(l)/(omega(idm)*sqrt(eps0)*spsh(l)))**2 write(6,*)dos(idm,jx),taui(l),rl(l),spsh(l) 10 continue dos(idm,jx)=-igf*2.d0*omega(idm)*sqrt(eps0)*dos(idm,jx)/(2.d0*pi) write(6,*)dos(idm,jx),igf * * Factor 2.*omega*sqrt(eps0) here corrects for the fact that we * are looking and $\rho(\sg)$ which is determined from the * relation $\rho(E)dE=\rho(\sg)d\sg$, where dE=2\sg d\sg and * \rho(E)=- Im G/\pi * * On the output, dos(idm,jx) contains the LDOS per mode of the field. rmg=dabs(omega(idm)- edw) WRITE(NOUT+1+jx,*) rmg, dos(idm,jx) if (id.eq.0) go to 100 rmg=log(dabs(omega(1)- edw)) rmg=rmg-log(dabs(omega(2)- edw)) rmg=-dabs(rmg) if (idm.eq.1) then ddos(jx)=log(dos(1,jx))-log(dos(2,jx)) else ddos(jx)=log(dos(2,jx))-log(dos(1,jx)) end if WRITE(NOUT+1+jx,*)'eta=', ddos(jx)/rmg 100 continue *---------------------------------------------------------------------- write(6,*)'omega=',omega(idm) 200 continue close(nout+1) close(nout+2) close(nout+3) close(nout+4) close(nout+5) *oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo end C (C) Copr. X/1999 Alexander Moroz