2
! Copyright (C) 1996-2016 The SIESTA group
3
! This file is distributed under the terms of the
4
! GNU General Public License: see COPYING in the top directory
5
! or http://www.gnu.org/copyleft/gpl.txt .
6
! See Docs/Contributors.txt for a list of contributors.
9
SUBROUTINE STSOFR( NA, NO, NUO, MAXNA, NSPIN,
10
. ISA, IPHORB, INDXUO, LASTO, XA, CELL,
11
. RPSI, IPSI, E, INDW, NWF, NUMWF, NK, K,
12
. ARMUNI, IUNITCD, RMAXO,
13
. EMIN, NE, DELTAE, ETA, XSTS, BFUNC )
14
C **********************************************************************
15
C Computes the local density of states as a function of energy
16
C at a given point in real space, for the simulation of STS spectra
17
C Coded by P. Ordejon. July 2004
18
C **********************************************************************
24
USE LISTSC_MODULE, ONLY: LISTSC
28
INTEGER, INTENT(IN) ::
29
. NA, NO, NUO, IUNITCD, NK,
30
. NSPIN, MAXNA, NWF(NK), NUMWF,
31
. ISA(NA), IPHORB(NO), INDXUO(NO), LASTO(0:NA),
32
. INDW(NK,NUMWF), NE, BFUNC
33
real(dp), INTENT(IN) ::
34
. XA(3,NA), ARMUNI, RMAXO,
35
. EMIN, DELTAE, ETA, XSTS(3)
37
real(dp), INTENT(IN) ::
39
. RPSI(NUO,NK,NUMWF,NSPIN), IPSI(NUO,NK,NUMWF,NSPIN),
40
. E(NK,NUMWF,NSPIN), K(NK,3)
42
C ****** INPUT *********************************************************
43
C INTEGER NA : Total number of atoms in Supercell
44
C INTEGER NO : Total number of orbitals in Supercell
45
C INTEGER NUO : Total number of orbitals in Unit Cell
46
C INTEGER MAXNA : Maximum number of neighbours of any atom
47
C INTEGER NSPIN : Number of different spin polarizations
48
C Nspin = 1 => unpolarized, Nspin = 2 => polarized
49
C INTEGER ISA(NA) : Species index of each atom
50
C INTEGER IPHORB(NO) : Orital index of each orbital in its atom
51
C INTEGER INDXUO(NO) : Equivalent orbital in unit cell
52
C INTEGER LASTO(0:NA) : Last orbital of each atom in array iphorb
53
C REAL*8 XA(3,NA) : Atomic positions in cartesian coordinates
55
C REAL*8 CELL(3,3) : Supercell vectors CELL(IXYZ,IVECT)
57
C REAL*8 RPSI(NUO,NK,NUMWF,NSPIN): Wave function coefficients (real part)
58
C REAL*8 IPSI(NUO,NK,NUMWF,NSPIN): Wave function coefficients (imag part)
59
C REAL*8 E(NK,NUMWF,NSPIN) : Eigen energies
60
C INTEGER INDW(NUMWF,NK) : Index of the wavefunctions
61
C INTEGER NWF(NK) : Number of wavefncts to print for each k-point
62
C INTEGER NUMWF : Max num of wavefncts to print a given k-point
63
C INTEGER NK : Number of k-points
64
C REAL*8 K(NK,3) : k-points
65
C REAL*8 ARMUNI : Conversion factor for the charge density
66
C INTEGER IUNITCD : Unit of the charge density
67
C REAL*8 RMAXO : Maximum range of basis orbitals
68
C REAL*8 EMIN : Minimum energy in STS curve
69
C INTEGER NE : Number of points in STS curve
70
C REAL*8 DELTAE : Energy step between points in STS curve
71
C REAL*8 ETA : Lorenzian or gaussian width
72
C REAL*8 XSTS(3) : Position where STS is calculated
73
C INTEGER BFUNC : Broadening function
74
C 1 = Lorentzian 2 = Gaussian
75
C **********************************************************************
77
INTEGER, DIMENSION(:), ALLOCATABLE :: JNA
79
REAL, DIMENSION(:,:), allocatable :: RWF, IMWF, MWF, PWF
81
real(dp), DIMENSION(:), ALLOCATABLE :: R2IJ, LDOS
83
real(dp), DIMENSION(:,:), ALLOCATABLE :: XIJ
86
. IA, ISEL, NNA, UNIT, IE
89
. I, J, IN, IAT1, IAT2, JO, IO, IUO, IAVEC, IAVEC1, IAVEC2,
90
. IS1, IS2, IPHI1, IPHI2, IND, IX, IY, IZ, NX, NY, NZ, IWF,
94
. RMAX, XPO(3), RMAX2, XVEC1(3),
95
. XVEC2(3), PHIMU, GRPHIMU(3),
96
. OCELL(3,3), PHASE, SI, CO, PI, ENER, PSI2
99
. RWAVE, RWAVEUP, RWAVEDN,
100
. IWAVE, IWAVEUP, IWAVEDN
105
. SNAME*40, FNAME*60,
106
. CHAR1*10, CHAR2*10,
110
. IO_ASSIGN, IO_CLOSE, NEIGHB
112
C **********************************************************************
113
C INTEGER IA : Atom whose neighbours are needed.
114
C A routine initialization must be done by
115
C a first call with IA = 0
116
C If IA0=0, point X0 is used as origin instead
117
C INTEGER ISEL : Single-counting switch (0=No, 1=Yes). If ISEL=1,
118
C only neighbours with JA.LE.IA are included in JNA
119
C INTEGER NNA : Number of non-zero orbitals at a point in
121
C INTEGER JNA(MAXNA) : Atom index of neighbours. The neighbours
122
C atoms might be in the supercell
123
C REAL*8 XIJ(3,MAXNA) : Vectors from point in real space to orbitals
124
C REAL*8 R2IJ(MAXNA) : Squared distance to atomic orbitals
125
C REAL*8 XPO(3) : Coordinates of the point of the plane respect
126
C we are going to calculate the neighbours orbitals
127
C INTEGER IZA(NA) : Atomic number of each atom
128
C **********************************************************************
131
C Allocate some variables ---------------------------------------------
132
IF (NSPIN .NE. 1 .AND. NSPIN .NE. 2) THEN
133
WRITE(6,*)'BAD NUMBER NSPIN IN STS.F'
134
WRITE(6,*)'NSPIN = ',NSPIN
135
WRITE(6,*)'IT MUST BE 1 => NON POLARIZED, OR 2 = > POLARIZED'
139
PI = 4.0D0 * ATAN(1.0D0)
142
CALL MEMORY('A','D',SIZE(LDOS),'stsofr')
143
C Initiallize ldos vector ---------------------------------------------
149
C Initialize neighbour subroutine --------------------------------------
154
IF (ALLOCATED(JNA)) THEN
155
CALL MEMORY('D','I',SIZE(JNA),'stsofr')
158
IF (ALLOCATED(R2IJ)) THEN
159
CALL MEMORY('D','D',SIZE(R2IJ),'stsofr')
162
IF (ALLOCATED(XIJ)) THEN
163
CALL MEMORY('D','D',SIZE(XIJ),'stsofr')
167
CALL MEMORY('A','I',MAXNA,'stsofr')
168
ALLOCATE(R2IJ(MAXNA))
169
CALL MEMORY('A','D',MAXNA,'stsofr')
170
ALLOCATE(XIJ(3,MAXNA))
171
CALL MEMORY('A','D',3*MAXNA,'stsofr')
178
CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL,
179
. NNA, JNA, XIJ, R2IJ, FIRST )
185
C Assign the point where STS will be computed
191
C Open file to store sts spectrum -----------------------------------
192
SNAME = FDF_STRING('SystemLabel','siesta')
193
FNAME = TRIM(SNAME)//'.STS'
195
OPEN(UNIT = UNIT, FILE = FNAME, STATUS = 'UNKNOWN',
196
. FORM = 'FORMATTED')
201
. ' Generating STS spectra at point'
203
. XPO(1),XPO(2),XPO(3)
205
. ' Please, wait....'
208
C Localize non-zero orbitals each point in real space ---------------
214
CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL,
215
. NNA, JNA, XIJ, R2IJ, FIRST )
218
C Loop over wavefunctions
226
C Initialize the value of the wave function -----------------------
234
C Loop over Non-zero orbitals ------------------------------------------
236
IF( R2IJ(IAT1) .GT. RMAX2 ) EXIT
240
XVEC1(1) = -XIJ(1,IAT1)
241
XVEC1(2) = -XIJ(2,IAT1)
242
XVEC1(3) = -XIJ(3,IAT1)
244
C Calculate phase of Bloch wavefunctions
246
PHASE = K(IK,1)*(XPO(1)+XIJ(1,IAT1))+
247
. K(IK,2)*(XPO(2)+XIJ(2,IAT1))+
248
. K(IK,3)*(XPO(3)+XIJ(3,IAT1))
253
DO 120 IO = LASTO(IAVEC1-1) + 1, LASTO(IAVEC1)
256
CALL PHIATM( IS1, IPHI1, XVEC1, PHIMU, GRPHIMU )
258
IF ( NSPIN .EQ. 1 ) THEN
259
RWAVE = RWAVE + PHIMU *
260
. (RPSI(IUO,IK,IWF,1)*CO - IPSI(IUO,IK,IWF,1)*SI)
261
IWAVE = IWAVE + PHIMU *
262
. (RPSI(IUO,IK,IWF,1)*SI + IPSI(IUO,IK,IWF,1)*CO)
263
ELSEIF (NSPIN .EQ. 2) THEN
264
RWAVEUP = RWAVEUP + PHIMU *
265
. (RPSI(IUO,IK,IWF,1)*CO - IPSI(IUO,IK,IWF,1)*SI)
266
IWAVEUP = IWAVEUP + PHIMU *
267
. (RPSI(IUO,IK,IWF,1)*SI + IPSI(IUO,IK,IWF,1)*CO)
268
RWAVEDN = RWAVEDN + PHIMU *
269
. (RPSI(IUO,IK,IWF,2)*CO - IPSI(IUO,IK,IWF,2)*SI)
270
IWAVEDN = IWAVEDN + PHIMU *
271
. (RPSI(IUO,IK,IWF,2)*SI + IPSI(IUO,IK,IWF,2)*CO)
280
C Loop over energies to calculate the local density of states
286
IF (NSPIN .EQ. 1) THEN
287
PSI2 = RWAVE**2+IWAVE**2
288
IF (BFUNC .EQ. 1) THEN
289
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
290
. (ETA/PI) / ((ENER - E(IK,IWF,1))**2 + ETA**2)
291
ELSE IF (BFUNC .EQ. 2) THEN
292
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
293
. (1.0d0/(ETA*DSQRT(PI))) *
294
. DEXP(-((ENER - E(IK,IWF,1))/ETA)**2)
296
STOP 'Wrong convolution function in stsofr.f'
298
ELSE IF (NSPIN .EQ. 2) THEN
299
PSI2 = RWAVEUP**2+IWAVEUP**2
300
IF (BFUNC .EQ. 1) THEN
301
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
302
. (ETA/PI) / ((ENER - E(IK,IWF,1))**2 + ETA**2)
303
ELSE IF (BFUNC .EQ. 2) THEN
304
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
305
. (1.0d0/(ETA*DSQRT(PI))) *
306
. DEXP(-((ENER - E(IK,IWF,1))/ETA)**2)
308
STOP 'Wrong convolution function in stsofr.f'
310
PSI2 = RWAVEDN**2+IWAVEDN**2
311
IF (BFUNC .EQ. 1) THEN
312
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
313
. (ETA/PI) / ((ENER - E(IK,IWF,2))**2 + ETA**2)
314
ELSE IF (BFUNC .EQ. 2) THEN
315
LDOS(IE) = LDOS(IE) + (ARMUNI/NK) * 2.0 * PSI2 *
316
. (1.0d0/(ETA*DSQRT(PI))) *
317
. DEXP(-((ENER - E(IK,IWF,2))/ETA)**2)
319
STOP 'Wrong convolution function in stsofr.f'
332
C Print out LDOS -------------------------------------------------
333
WRITE(UNIT,'(A)') '# LDOS vs. Energy curve'
334
WRITE(UNIT,'(A)') '# '
335
WRITE(UNIT,'(A,3e13.5,A)')
336
. '# LDOS computed on point ',(XSTS(I),I=1,3),' bohr'
337
IF (IUNITCD .EQ. 1) THEN
338
WRITE(UNIT,'(A)') '# in units of elecron/(bohr)**3/eV'
339
ELSE IF (IUNITCD .EQ. 2) THEN
340
WRITE(UNIT,'(A)') '# in units of elecron/(ang)**3/eV'
341
ELSE IF (IUNITCD .EQ. 3) THEN
342
WRITE(UNIT,'(A)') '# in units of elecron/unitcell/eV'
344
STOP 'Wrong IUNITCD in stsofr.f'
346
WRITE(UNIT,'(A)') '# '
347
WRITE(UNIT,'(A)') '# Energy in eV)'
348
WRITE(UNIT,'(A)') '# '
349
WRITE(UNIT,'(A)') '# Energy LDOS'
350
WRITE(UNIT,'(A)') '# '
354
WRITE(UNIT,'(2f15.5)') ENER,LDOS(IE)
360
. ' Your STS output file is:'
361
WRITE(6,'(A,A)') ' ',FNAME
366
CALL MEMORY('D','D',SIZE(LDOS),'stsofr')
added
removed
Util/Denchar/Src/stsofr.f
