« back to all changes in this revision
Viewing changes to Util/STM/simple-stm/plstm.f90
- browse files at revision 605
- compare with another revision
- download diff
- download tarball
- view history from revision 605
- Committer: Alberto Garcia
- Date: 2019-09-02 14:09:43 UTC
- mfrom: (427.6.323 trunk)
- Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)
- files added:
- Src/bloch_unfold.F90
- Tests/bi2se3_3ql
- Util/COOP/Tests/bi2se3_3ql
- Util/COOP/Tests/ge_fatbands_so
- Util/STM/ol-stm/Examples/Benzene
- Util/STM/ol-stm/Examples/ML-Cr
- files removed:
- README_TRANSIESTA
- files renamed:
- Util/STM/simple-stm/plstm.f => Util/STM/simple-stm/plstm.f90
- Util/WFS/readwfx.f => Util/WFS/readwfx.f90
- files modified:
- Docs/Contributors.txt
added
removed
Util/STM/simple-stm/plstm.f90
1
1
!
2
! Copyright (C) 1996-2016 The SIESTA group
2
! Copyright (C) 1996-2016 The SIESTA group
3
3
! This file is distributed under the terms of the
4
4
! GNU General Public License: see COPYING in the top directory
5
5
! or http://www.gnu.org/copyleft/gpl.txt.
6
6
! See Docs/Contributors.txt for a list of contributors.
7
7
!
8
8
9
program plstm
10
11
c****************************************************************************
12
c PLSTM version 1.1.2
13
c
14
c This program PLSTM reads a charge density or local density of states
15
c generated by SIESTA, and simulates STM images at constant current
16
c or constant height mode using the Tersoff-Hamann approximation.
17
c
18
c Written by P. Ordejon, June 2001
19
c Last modification: June 2003
20
c (basic structure from Plrho package of J.M.Soler)
21
c****************************************************************************
22
c CAVEATS:
23
c
24
c This version works assuming that the scanning plane is the XY plane.
25
c This plane must be perpendicular to the third lattice vector of
26
c the supercell (the Z direction).
27
c****************************************************************************
28
c USAGE:
29
c
30
c This program reads files generated from SIESTA, with information of
31
c the charge density (filename.RHO) or local density of states (filename.LDOS)
32
c and computes a simulations of STM images, in the Tersoff-Hamann
33
c approximation. Two modes are available: constant current (simulated
34
c by computing a constant density surface z=z(x,y)) and constant hight
35
c (obtaining the charge density at the tip position at a given height).
36
c
37
c The program needs two input files:
38
c
39
c 1) Main input file, read by standard input. A sample of input file is:
40
c
41
c --- begin input file ---
42
c graf
43
c rho
44
c constant-current
45
c 1.d-4
46
c unformatted
47
c --- end input file ---
48
c
49
c where:
50
c - The first line is the label of the system, as in SIESTA SystemLabel.
51
c Files will be searched as SystemLabel.* (in the example, graf.RHO).
52
c - The second line is the task, which should be either rho or ldos
53
c (in lowercase!!).
54
c - The third line specifies the STM mode: 'constant-current' or
55
c 'constant-height' (in lowercase).
56
c - The fourth line is a value that determines the details of the
57
c STM image. In the case of 'constant-current' mode, this is the
58
c value of the density (in units of e/bohr**3). at which the isosurface
59
c is computed. For the 'constant-height' mode, the value specifies
60
c the Z (in bohr) level at which the charge is computed.
61
c - The fifth line indicates if the grid data file SystemLabel.TASK
62
c is formatted or unformatted (the latter being the standard option
63
c in SIESTA)
64
c
65
c 2) SystemLabel.TASK file: this is a file generated by SIESTA, with
66
c the values of the appropriate quantity on the grid.
67
c In example above: grid.RHO. You should copy it from the directory
68
c with your SIESTA output files.
69
c
70
c The program generates some informative output on standard output,
71
c and writes one file. In the case of 'constant-current' mode,
72
c this file's name is SystemLabel.CC.STM, and contains the X,Y,Z
73
c values of the isosurface (a grid of X,Y and the value of Z(X,Y) of
74
c the isosurface). In the case of 'constant-height' mode, the name
75
c is SystemLabel.CH.STM, and contains the values X,Y,RHO for each
76
c X,Y of the grid, where RHO is the charge computed at the point X,Y,Z
77
c (Z being the hight specified in the input).
78
c
79
c If your SIESTA calculation was spin polarized, the program adds
80
c the two spin components, so the output represents the total charge
81
c density.
82
c****************************************************************************
83
84
85
86
implicit none
87
88
integer, parameter :: dp = selected_real_kind(10,100)
89
c Internal parameters
90
c maxp : Maximun number of points
91
integer, parameter :: maxp = 10000000
92
93
c Internal variables
94
character
95
. name*75, fform*12, fname*80, oname*80, task*15,
96
. mode*25
97
integer
98
. i, ip, is, j, mesh(3), np, nspin, nt, Ind, iz, iy
99
100
real :: f(maxp,2), fvalue, rho
101
real(dp) :: cell(3,3)
102
103
c Read plot data
104
read(5,*) name
105
read(5,*) task
106
read(5,*) mode
107
read(5,*) fvalue
108
read(5,*) fform
109
110
c Read density
111
112
if (task .eq. 'ldos') then
113
fname = trim(name)//'.LDOS'
114
else if (task .eq. 'rho') then
115
fname = trim(name)//'.RHO'
116
else
117
stop ' ERROR: Task should be RHO or LDOS'
9
program plstm
10
11
12
! Program PLSTM reads a charge density or local density of states
13
! generated by SIESTA, and simulates STM images at constant-current
14
! or constant-height mode using the Tersoff-Hamann approximation.
15
16
! Written by P. Ordejon, June 2001
17
! (basic structure from Plrho package of J.M.Soler)
18
! Spin capability and restructuring by Alberto García (March 2019)
19
20
! CAVEATS:
21
!
22
! This version works assuming that the scanning plane is the XY plane.
23
! This plane must be perpendicular to the third lattice vector of
24
! the supercell (the Z direction).
25
26
! USAGE:
27
28
! This program reads grid files generated from SIESTA or other programs
29
! (such as ol-stm/wfs2ldos), with information on the local density of
30
! states (filename.LDOS) and computes a simulated STM image, in the
31
! Tersoff-Hamann approximation. Two modes are available:
32
! constant-current (simulated by computing a constant ldos surface
33
! z=z(x,y)) and constant-height (obtaining the ldos at the tip position
34
! at a given height).
35
36
! The program uses command-line options. Type 'plstm -h' for usage.
37
38
! It is possible to analyze just the 'charge', or arbitrary components
39
! of the spin, by using the '-s' or '-v' options.
40
! '-v', followed by a line with three real numbers,
41
! serves to simulate a spin-polarized tip. The scalar product of the
42
! "tip spin vector" and the "local spin vector" is recorded at each point.
43
44
! The program generates some informative output on standard output,
45
! including valid ranges, and writes one file. In the case of
46
! 'constant-current' mode, this file's name is SystemLabel.*.CC.STM,
47
! and contains the X,Y,Z values of the isosurface (a grid of X,Y and
48
! the value of Z(X,Y) of the isosurface). In the case of
49
! 'constant-height' mode, the name is SystemLabel.*.CH.STM, and
50
! contains the values X,Y,RHO for each X,Y of the grid, where RHO is
51
! the charge (or spin component) computed at the point X,Y,Z (Z being
52
! the height specified in the input).
53
54
! In the filenames, '*' stands for the 'spin code'
55
! (as entered with the '-s' flag).
56
57
use m_gridfunc, only: monoclinic_z
58
use m_gridfunc, only: gridfunc_t, read_gridfunc
59
use m_getopts
60
61
implicit none
62
63
type(gridfunc_t) :: gf
64
65
integer, parameter :: dp = selected_real_kind(10,100)
66
67
character(len=200) :: opt_arg
68
character(len=10) :: opt_name
69
integer :: nargs, iostat, n_opts, nlabels
70
71
logical :: debug = .false.
72
73
real, allocatable :: rho(:), f2d(:,:)
74
75
character name*75, fform*12, fname*80, oname*80, task*15, &
76
mode*25, spin_code*1
77
integer :: i, ip, is, j, mesh(3), np, nspin, nt, Ind, iz, iy
78
integer :: k1, k2, ic, ix
79
real :: fvalue, sx, sy, sz, x, y
80
real :: fmin, fmax
81
real(dp) :: zmin, zmax, stepz
82
integer :: n3, q1, q2
83
real(dp) :: cell(3,3), origin(3)
84
real(dp) :: dxdm(2,2) ! 2D plane grid vectors
85
real(dp) :: tip_spin(3)
86
logical :: get_current_range, get_height_range
87
integer :: nx, ny
88
!
89
! Process options
90
!
91
! defaults
92
get_current_range = .false.
93
get_height_range = .false.
94
spin_code = 'q' ! default value
95
nx = 1
96
ny = 1
97
n_opts = 0
98
do
99
call getopts('hdz:i:s:v:X:Y:IH',opt_name,opt_arg,n_opts,iostat)
100
if (iostat /= 0) exit
101
select case(opt_name)
102
case ('d')
103
debug = .true.
104
case ('z')
105
mode = 'constant-height'
106
read(opt_arg,*) fvalue
107
case ('i')
108
mode = 'constant-current'
109
read(opt_arg,*) fvalue
110
case ('s')
111
read(opt_arg,*) spin_code
112
case ('v')
113
spin_code = 'v'
114
read(opt_arg,*) tip_spin(1:3)
115
case ('X')
116
read(opt_arg,*) nx
117
case ('Y')
118
read(opt_arg,*) ny
119
case ('I')
120
get_current_range = .true.
121
case ('H')
122
get_height_range = .true.
123
case ('h')
124
call manual()
125
STOP
126
case ('?',':')
127
write(0,*) "Invalid option: ", opt_arg(1:1)
128
write(0,*) "Use -h option for manual"
129
write(0,*) ""
130
call manual()
131
STOP
132
end select
133
enddo
134
135
nargs = command_argument_count()
136
nlabels = nargs - n_opts + 1
137
if (nlabels /= 1) then
138
write(0,*) "Use -h option for manual"
139
write(0,*) ""
140
call manual()
141
STOP
142
endif
143
144
call get_command_argument(n_opts,value=fname,status=iostat)
145
if ( iostat /= 0 ) then
146
stop "Cannot get LDOS file"
147
end if
148
149
if (spin_code == 'q') then
150
write(0,*) "Using 'total charge' ('q') mode"
151
endif
152
153
write(6,*) 'Reading grid data from file ',trim(fname)
154
155
!! SHOULD MAKE SURE OF THE UNITS... although for plots
156
! it should not make much difference
157
158
! We might want to try to detect file extension and use
159
! the proper reader.
160
call read_gridfunc(fname,gf)
161
162
np = product(gf%n(1:3))
163
nspin = gf%nspin
164
cell = gf%cell
165
mesh = gf%n
166
origin = gf%origin
167
168
write(6,*)
169
write(6,*) 'Cell vectors (bohr)'
170
write(6,*)
171
write(6,*) cell(:,1)
172
write(6,*) cell(:,2)
173
write(6,*) cell(:,3)
174
write(6,*)
175
write(6,*) 'Grid mesh: ',mesh(1),'x',mesh(2),'x',mesh(3)
176
write(6,*)
177
write(6,*) 'nspin = ',nspin
178
write(6,*)
179
write(6,"(a,3f10.5)") 'Box origin (bohr): ', origin(1:3)
180
181
if (.not. monoclinic_z(cell)) then
182
write(0,*) 'The cell is not monoclinic with ' // &
183
'c lattice vector along z...'
184
stop " *** Unsuitable cell for STM program"
118
185
endif
119
186
120
write(6,*)
121
write(6,*) '**********************************************'
122
if (mode .eq. 'constant-current') then
123
write(6,*) 'Calculating STM image in Constant Current mode'
124
write(6,*) 'The STM image is obtained as the isosurface of'
125
write(6,*) 'constant charge density RHO =', fvalue,' e/Bohr**3'
126
oname = trim(name)//'.CC.STM'
127
else if (mode .eq. 'constant-height') then
128
write(6,*) 'Calculating STM image in Constant Height mode'
129
write(6,*) 'The STM image is obtained as the value of the'
130
write(6,*) 'charge at a given tip height Z = ', fvalue, 'Bohr'
131
oname = trim(name)//'.CH.STM'
187
! Note that the last slice with information in the file is
188
! the (n3-1)th plane
189
n3 = mesh(3)
190
zmax = (n3-1) * cell(3,3) / n3 + origin(3)
191
zmin = origin(3)
192
if (n3==1) then
193
write(6,"(/,a)") "The file contains a single plane of data"
194
write(6,"(/,a)") "Only constant-height mode is allowed"
195
stepz = 0.0
196
else
197
stepz = cell(3,3) / n3
198
endif
199
write(6,"(a,3f12.5)") "Zmin, Zmax (bohr): ", zmin, zmax
200
201
202
if (mode .eq. 'constant-current') then
203
204
if (n3==1) then
205
STOP "Constant-current mode not available for n3=1"
206
endif
207
208
! It does not make physical sense to get topography
209
! for spin components
210
if (spin_code /= "q") then
211
STOP "Constant-current mode only available " // &
212
"for 'q' (charge)"
213
endif
214
215
else if (mode .eq. 'constant-height') then
216
217
if (n3==1) then
218
write(6,"(/,a)") "The file contains a single plane of data"
219
write(6,"(/,a)") "Height set to plane's height"
220
fvalue = zmin
221
else
222
223
! Check that value is within range
224
if ((fvalue > zmax) .or. (fvalue < zmin)) then
225
STOP 'Z requested is beyond box limits'
226
endif
227
endif
228
132
229
else
133
write(6,*) 'ERROR: mode must be either constant current'
134
write(6,*) ' or constant height (in lower case)'
230
write(6,*) 'ERROR: mode must be either constant-current'
231
write(6,*) ' or constant-height (in lower case)'
135
232
stop
136
233
endif
137
write(6,*) '**********************************************'
138
write(6,*)
139
140
141
write(6,*)
142
write(6,*) 'Reading grid data from file ',trim(fname)
143
144
open( unit=1, file=fname, status='old', form=fform )
145
if (fform .eq. 'formatted') then
146
read(1,*) cell
147
read(1,*) mesh, nspin
148
np = mesh(1) * mesh(2) * mesh(3)
149
if (np .gt. maxp) stop 'plrho: parameter maxp too small'
150
do is=1,nspin
151
Ind=0
152
do iz=1,mesh(3)
153
do iy=1,mesh(2)
154
read(1,'(e15.6)') (f(Ind+ip,is),ip=1,mesh(1))
155
Ind=Ind+mesh(1)
156
enddo
157
enddo
158
enddo
159
else
160
read(1) cell
161
read(1) mesh, nspin
162
np = mesh(1) * mesh(2) * mesh(3)
163
C write(6,*) cell,mesh,nspin
164
if (np .gt. maxp) stop 'plrho: parameter maxp too small'
165
do is=1,nspin
166
Ind=0
167
do iz=1,mesh(3)
168
do iy=1,mesh(2)
169
read(1) (f(Ind+ip,is),ip=1,mesh(1))
170
Ind=Ind+mesh(1)
171
enddo
172
enddo
173
enddo
174
endif
175
close(1)
176
177
write(6,*)
178
write(6,*) 'Cell vectors'
179
write(6,*)
180
write(6,*) cell(1,1),cell(2,1),cell(3,1)
181
write(6,*) cell(1,2),cell(2,2),cell(3,2)
182
write(6,*) cell(1,3),cell(2,3),cell(3,3)
183
write(6,*)
184
write(6,*) 'Grid mesh: ',mesh(1),'x',mesh(2),'x',mesh(3)
185
write(6,*)
186
write(6,*) 'nspin = ',nspin
187
write(6,*)
188
189
190
c Get total density if spin polarized
191
if (nspin .eq. 2) then
192
do ip = 1,np
193
rho = f(ip,1) + f(ip,2)
194
f(ip,2) = f(ip,1) - f(ip,2)
195
f(ip,1) = rho
196
enddo
197
endif
198
199
c Generate x,y,z surface (dump to output file)
200
201
if (mode .eq. 'constant-current') then
202
call isocharge( cell, mesh, mesh, f(1,1), fvalue, oname)
203
else
204
call isoz( cell, mesh, mesh, f(1,1), fvalue, oname)
205
endif
206
207
end
208
209
210
SUBROUTINE ISOCHARGE( CELL, NMESH, NSPAN, F, FVALUE, ONAME)
211
212
C *******************************************************************
213
C Calculates the surface z=z(x,y) with constant function value.
214
C The surface is determined by the condition function=value, and
215
C it is printed in a file as x,y,z. The function must
216
C be given in a regular 3-D grid of points.
217
C Notice single precision in this version
218
C
219
C Written by P. Ordejon. June 2001.
220
C from plsurf.f (written by J. M. Soler)
221
C ************************* INPUT ***********************************
222
C REAL*8 CELL(3,3) : Unit cell vectors CELL(ixyz,ivector)
223
C INTEGER NMESH(3) : Number of mesh divisions of each vector
224
C INTEGER NSPAN(3) : Physical dimensions of array F
225
C See usage section for more information
226
C REAL F(*) : Function such that F=FVALUE determines
227
C the shape of the solid surface.
228
C REAL FVALUE : Value such that F=FVALUE
229
C determines the shape of the solid surface.
230
C CHARACTER*80 ONAME : Output file name
231
C ************************* OUTPUT **********************************
232
C None. Results are printed on ONAME file
233
C *******************************************************************
234
235
236
C Next line is nonstandard but may be suppressed
234
write(6,*)
235
236
allocate(rho(product(mesh(1:3))))
237
238
call get_function(nspin, spin_code, mesh, gf%val, &
239
rho, tip_spin, fmin, fmax)
240
241
select case (spin_code)
242
case ( 'q' )
243
write(6,*) "Using 'total charge' ('q') mode"
244
case ( 'x' )
245
write(6,*) "Using 'x-component of spin' ('x') mode"
246
case ( 'y' )
247
write(6,*) "Using 'y-component of spin' ('y') mode"
248
case ( 'z' )
249
write(6,*) "Using 'z-component of spin' ('z') mode"
250
case ( 's' )
251
write(6,*) "Using 'total spin' ('s') mode"
252
case ( 'v' )
253
write(6,"(a,3f10.5)") "Using a 'polarized tip' ('-v') " // &
254
"with spin: ", tip_spin(1:3)
255
end select
256
257
write(6,"(a,3g20.8,/)") "Range of values of processed function: ", fmin, fmax
258
259
allocate(f2d(0:mesh(1)-1,0:mesh(2)-1))
260
261
if (mode .eq. 'constant-current') then
262
write(6,"(/a)") 'Calculating STM image in Constant Current mode'
263
write(6,"(a)") 'The STM image is obtained as the isosurface of'
264
write(6,"(a)") 'constant function value =', fvalue,' e/Bohr**3'
265
! Check that value is within range
266
if ((fvalue > fmax) .or. (fvalue < fmin)) then
267
STOP 'Function value outside range of values in box'
268
endif
269
oname = trim(fname)// "." // spin_code // '.CC.STM'
270
271
! Generate x,y,z surface (dump to output file)
272
call isocharge( stepz, origin, mesh, rho, fvalue, f2d)
273
274
else
275
276
! Dump slice to output file
277
write(6,"(/a)") 'Calculating STM image in Constant Height mode'
278
write(6,"(a)") 'The STM image is obtained as the value of the'
279
write(6,"(a,f9.4,a)") 'charge at a given tip height Z = ', fvalue, 'Bohr'
280
oname = trim(fname)// "." // spin_code // '.CH.STM'
281
282
call isoz( stepz, origin, mesh, rho, fvalue, f2d )
283
endif
284
285
! Write 2D file info
286
287
OPEN( unit=2, file=oname )
288
write(6,*) 'Writing STM image in file ', trim(oname)
289
290
DO IC = 1,2
291
DO IX = 1,2
292
DXDM(IX,IC) = CELL(IX,IC) / MESH(IC)
293
ENDDO
294
ENDDO
295
296
! Possibly use nx, ny multipliers here
297
do k1 = 0, nx*mesh(1) - 1
298
do k2 = 0, ny*mesh(2) - 1
299
x = dxdm(1,1) * k1 + dxdm(1,2) * k2
300
y = dxdm(2,1) * k1 + dxdm(2,2) * k2
301
q1 = mod(k1,mesh(1))
302
q2 = mod(k2,mesh(2))
303
write(2,*) x, y, f2d(q1,q2)
304
enddo
305
write(2,*)
306
enddo
307
close(2)
308
deallocate(rho,f2d)
309
310
CONTAINS
311
subroutine manual()
312
write(0,"(a)") " -------------------"
313
write(0,"(a)") " Usage: plstm [options] LDOSfile"
314
write(0,"(a)") " "
315
write(0,"(a)") " "
316
write(0,"(a)") " OPTIONS: "
317
write(0,"(a)") " "
318
write(0,"(a)") " -h Print this help"
319
write(0,"(a)") " -d Print debugging info"
320
write(0,"(a)") " "
321
write(0,"(a)") " -i current Constant-current calculation"
322
write(0,"(a)") " with 'current' in e/bohr**3 "
323
write(0,"(a)") " -z height Constant-height calculation"
324
write(0,"(a)") " with 'height' in bohr "
325
write(0,"(a)") " "
326
write(0,"(a)") " -s {q,x,y,z,s} Spin code (default 'q' for total 'charge')"
327
write(0,"(a)") " (x|y|z) select cartesian components of spin "
328
write(0,"(a)") " s selects total spin magnitude "
329
write(0,"(a)") " -v 'ux uy uz' Tip spin direction for selection of spin component"
330
write(0,"(a)") " The vector (ux,uy,uz) should be normalized"
331
write(0,"(a)") " "
332
write(0,"(a)") " -X NX Request multiple copies of plot domain along X"
333
write(0,"(a)") " -Y NY Request multiple copies of plot domain along Y"
334
write(0,"(a)") " "
335
write(0,"(a)") " -H Return range of height (not implemented yet)"
336
write(0,"(a)") " -I Return range of current (not implemented yet)"
337
write(0,"(a)") " -------------------"
338
339
end subroutine manual
340
end program plstm
341
342
subroutine get_function(nspin, spin_code, mesh, f, rho, &
343
tip_spin, fmin, fmax)
344
345
integer, parameter :: dp = selected_real_kind(10,100)
346
integer, parameter :: sp = kind(1.0)
347
348
integer, intent(in) :: nspin, mesh(3)
349
real(sp), intent(in) :: f(mesh(1)*mesh(2)*mesh(3),nspin)
350
character(len=1), intent(in) :: spin_code
351
real(dp), intent(in) :: tip_spin(3)
352
353
real(sp), intent(out) :: rho(mesh(1)*mesh(2)*mesh(3))
354
real(sp), intent(out) :: fmin, fmax
355
356
real :: sx, sy, sz
357
integer :: i, np
358
359
np = product(mesh(1:3))
360
361
if (nspin == 1) then
362
select case (spin_code)
363
case ( 'x', 'y', 'z', 's')
364
stop "Cannot choose spin for spin-less file"
365
case ( 'v')
366
stop "Cannot use tip spin for spin-less file"
367
case ( 'q')
368
rho(:) = f(:,1)
369
end select
370
else if (nspin == 2) then
371
select case (spin_code)
372
case ( 'x', 'y')
373
stop "Cannot choose x, y comps for collinear file"
374
case ( 'v')
375
stop "Cannot use tip spin for collinear file"
376
case ( 'q')
377
rho(:) = f(:,1) + f(:,2)
378
case ( 'z')
379
rho(:) = f(:,1) - f(:,2)
380
case ( 's')
381
rho(:) = abs(f(:,1) - f(:,2))
382
end select
383
else if (nspin == 4) then
384
select case (spin_code)
385
case ( 'x' )
386
rho(:) = 2.0 * f(:,3)
387
case ( 'y' )
388
rho(:) = 2.0 * f(:,4)
389
case ( 'z')
390
rho(:) = f(:,1) - f(:,2)
391
case ( 'q')
392
rho(:) = f(:,1) + f(:,2)
393
case ( 's')
394
do i = 1, np
395
sx = 2.0 * f(i,3)
396
sy = 2.0 * f(i,4)
397
sz = f(i,1) - f(i,2)
398
rho(i) = sqrt(sx**2 + sy**2 + sz**2)
399
enddo
400
case ( 'v')
401
do i = 1, np
402
sx = 2.0 * f(i,3)
403
sy = 2.0 * f(i,4)
404
sz = f(i,1) - f(i,2)
405
rho(i) = tip_spin(1) * sx + &
406
tip_spin(2) * sy + &
407
tip_spin(3) * sz
408
enddo
409
end select
410
endif
411
fmin = minval(rho)
412
fmax = maxval(rho)
413
end subroutine get_function
414
415
! Some extra optimizations and clarifications are possible in these routines
416
! for a monoclinic cell
417
418
SUBROUTINE ISOCHARGE( stepz, origin, NMESH, F, FVALUE, f2d)
419
420
! *******************************************************************
421
! Calculates the surface z=z(x,y) with constant function value.
422
! The surface is determined by the condition function=value, and
423
! it is printed in a file as x,y,z. The function must
424
! be given in a regular 3-D grid of points.
425
! Notice single precision in this version
426
427
! Written by P. Ordejon. June 2001.
428
! from plsurf.f (written by J. M. Soler)
429
! Modified by A. Garcia, March 2019
430
! ************************* INPUT ***********************************
431
! REAL(DP) STEPZ
432
! INTEGER NMESH(3) : Number of mesh divisions of each vector
433
! REAL F(:,:,:) : Function such that F=FVALUE determines
434
! the shape of the solid surface.
435
! REAL FVALUE : Value such that F=FVALUE
436
! determines the shape of the solid surface.
437
! ************************* OUTPUT **********************************
438
! f2d function
439
! *******************************************************************
440
237
441
IMPLICIT NONE
238
239
442
integer, parameter :: dp = selected_real_kind(10,100)
240
241
C Argument types and dimensions
242
INTEGER
243
. NMESH(3), NSPAN(3), NT
244
REAL(dp) :: CELL(3,3)
245
REAL
246
. F(*), FVALUE
247
CHARACTER
248
. ONAME*80
249
250
251
C Local variables and arrays
252
LOGICAL
253
. HIGH, ZERO
254
INTEGER
255
. IC, IP, IPM, IX, K1, K2, K3
256
REAL
257
. DXDM(3,3), ZK3
258
259
260
OPEN( unit=2, file=oname)
261
262
write(6,*) 'Writing STM image in file', trim(oname)
263
264
265
266
C Find Jacobian matrix dx/dmesh and its inverse
267
DO IC = 1,3
268
DO IX = 1,3
269
DXDM(IX,IC) = CELL(IX,IC) / NMESH(IC)
270
ENDDO
271
ENDDO
272
443
444
INTEGER, intent(in) :: NMESH(3)
445
REAL(DP), intent(in) :: STEPZ
446
REAL(DP), intent(in) :: origin(3)
447
REAL, intent(in) :: F(0:nmesh(1)-1,0:nmesh(2)-1,0:nmesh(3)-1)
448
REAL, intent(in) :: FVALUE
449
REAL, intent(out) :: F2D(0:nmesh(1)-1,0:nmesh(2)-1)
450
451
! Local variables and arrays
452
LOGICAL HIGH, ZERO
453
INTEGER IC, IX, K1, K2, K3
454
REAL(dp) ZK3
455
REAL f_cur, f_up
273
456
274
457
ZERO = .FALSE.
275
C Loop on mesh points
458
276
459
DO K1 = 0,NMESH(1)-1
277
DO K2 = 0,NMESH(2)-1
278
C z-direction is scanned from top to bottom
279
DO K3 = NMESH(3)-1,0,-1
280
281
C Check if all cube vertices are above or below equi-surface.
282
HIGH = .FALSE.
283
284
C Find mesh index of this point
285
IP = 1 + K1 + NSPAN(1) * K2 + NSPAN(1) * NSPAN(2) * K3
286
287
C Find if this point is above FVALUE
288
HIGH = (F(IP) .GT. FVALUE)
289
290
if (HIGH) then
291
if (K3 .eq. NMESH(3)-1) stop 'Surface above box boundary!!!'
292
293
C Linear interpolation to find z-coordinate of surface
294
295
IPM = 1 + K1 + NSPAN(1)*K2 + NSPAN(1)*NSPAN(2)*(K3+1)
296
ZK3 = (K3+1) - (FVALUE - F(IPM)) / (F(IP) - F(IPM))
297
298
write(2,*)
299
. ( DXDM(IX,1) * K1 +
300
. DXDM(IX,2) * K2 +
301
. DXDM(IX,3) * ZK3 , IX=1,3)
302
303
goto 10
304
endif
305
306
ENDDO
307
write(2,*)
308
. ( DXDM(IX,1) * K1 +
309
. DXDM(IX,2) * K2 +
310
. DXDM(IX,3) * ZK3 , IX=1,3)
311
ZERO = .TRUE.
312
313
10 ENDDO
314
write(2,*) " "
315
ENDDO
316
317
CLOSE(2)
460
DO K2 = 0,NMESH(2)-1
461
! z-direction is scanned from top to bottom
462
! We assume that the function decreases as z increases
463
! This is reasonable for the charge density, but not
464
! necessarily true for spin components.
465
! We might need to work with absolute values
466
DO K3 = NMESH(3)-1,0,-1
467
468
!Find if this point is above FVALUE
469
f_cur = f(k1,k2,k3)
470
HIGH = ( f_cur .GT. FVALUE)
471
472
if (HIGH) then ! our point is between k3 and k3+1
473
if (K3 .eq. NMESH(3)-1) then
474
stop 'Surface above box boundary!!!'
475
endif
476
! Linear interpolation to find z-coordinate of surface
477
f_up = f(k1,k2,k3+1)
478
! This is a real number in [k3,k3+1]
479
ZK3 = (K3+1) - (FVALUE - f_up) / (f_cur - f_up)
480
f2d(k1,k2) = stepz * zk3 + origin(3)
481
goto 10
482
endif
483
484
ENDDO
485
! Note that K3 is zero here, at the end of the loop
486
f2d(k1,k2) = origin(3)
487
ZERO = .TRUE.
488
489
10 ENDDO
490
ENDDO
318
491
319
492
IF (ZERO) THEN
320
write(6,*) 'WARNING: I could not find the isosurface'
321
write(6,*) ' for some X,Y points. For these, I set Z = 0'
493
write(6,*) 'WARNING: I could not find the isosurface'
494
write(6,*) ' for some X,Y points. For these, Z = ZMIN (0)'
322
495
ENDIF
323
496
324
END
325
326
327
SUBROUTINE ISOZ( CELL, NMESH, NSPAN, F, ZVALUE, ONAME)
497
END SUBROUTINE ISOCHARGE
498
499
SUBROUTINE ISOZ( stepz, origin, NMESH, F, ZVALUE, f2d)
328
500
329
C *******************************************************************
330
C Calculates the value of the function F at the plane Z=ZVALUE
331
C The function must be given in a regular 3-D grid of points.
332
C Notice single precision in this version
333
C
334
C Written by P. Ordejon. June 2001.
335
C from plsurf.f (written by J. M. Soler)
336
C ************************* INPUT ***********************************
337
C REAL*8 CELL(3,3) : Unit cell vectors CELL(ixyz,ivector)
338
C INTEGER NMESH(3) : Number of mesh divisions of each vector
339
C INTEGER NSPAN(3) : Physical dimensions of array F
340
C See usage section for more information
341
C REAL F(*) : Function such that F=FVALUE determines
342
C the shape of the solid surface.
343
C REAL ZVALUE : Z level where the function is written
344
C CHARACTER*80 ONAME : Output file name
345
C ************************* OUTPUT **********************************
346
C None. Results are printed on ONAME file
347
C *******************************************************************
348
349
350
C Next line is nonstandard but may be suppressed
501
! Calculates the value of the function F at the plane Z=ZVALUE
502
! The function must be given in a regular 3-D grid of points.
503
! Notice single precision in this version
504
505
! Written by P. Ordejon. June 2001.
506
! from plsurf.f (written by J. M. Soler)
507
! ************************* INPUT ***********************************
508
! REAL(DP) CELL(3,3) : Unit cell vectors CELL(ixyz,ivector)
509
! INTEGER NMESH(3) : Number of mesh divisions of each vector
510
! REAL F(*) : Function such that F=FVALUE determines
511
! the shape of the solid surface.
512
! REAL ZVALUE : Z level where the function is written
513
! CHARACTER*80 ONAME : Output file name
514
! ************************* OUTPUT **********************************
515
! real f2d(:,:)
516
! *******************************************************************
517
351
518
IMPLICIT NONE
352
519
353
520
integer, parameter :: dp = selected_real_kind(10,100)
354
C Argument types and dimensions
355
INTEGER
356
. NMESH(3), NSPAN(3), NT
357
REAL(dp) :: CELL(3,3)
358
REAL
359
. F(*), ZVALUE, Z, FV, ZM
360
CHARACTER
361
. ONAME*80
362
363
364
C Local variables and arrays
365
INTEGER
366
. IC, IP, IPM, IX, K1, K2, K3
367
REAL
368
. DXDM(3,3), ZK3
369
370
371
OPEN( unit=2, file=oname )
372
373
write(6,*) 'Writing STM image in file', trim(oname)
374
375
376
C Find Jacobian matrix dx/dmesh and its inverse
377
DO IC = 1,3
378
DO IX = 1,3
379
DXDM(IX,IC) = CELL(IX,IC) / NMESH(IC)
380
ENDDO
381
ENDDO
382
383
384
C Loop on mesh points
385
DO K1 = 0,NMESH(1)-1
386
DO K2 = 0,NMESH(2)-1
387
C z-direction is scanned from top to bottom
388
DO K3 = NMESH(3)-1,0,-1
389
390
C Find mesh index of this point
391
IP = 1 + K1 + NSPAN(1) * K2 + NSPAN(1) * NSPAN(2) * K3
392
393
C Calculate Z coordinate of this point:
394
Z = DXDM(3,1) * K1 + DXDM(3,2) * K2 + DXDM(3,3) * K3
395
396
IF (Z .LT. ZVALUE) THEN
397
398
C Linear interpolation to find the value of F at ZVALUE
399
400
IPM = 1 + K1 + NSPAN(1)*K2 + NSPAN(1)*NSPAN(2)*(K3+1)
401
ZM = DXDM(3,1) * K1 + DXDM(3,2) * K2 + DXDM(3,3) * (K3+1)
402
403
FV = F(IPM) + (F(IP)-F(IPM)) * (ZVALUE - ZM) / (Z - ZM)
404
405
write(2,*)
406
. ( DXDM(IX,1) * K1 +
407
. DXDM(IX,2) * K2 +
408
. DXDM(IX,3) * K3 , IX=1,2), FV
521
522
INTEGER, intent(in) :: NMESH(3)
523
REAL(dp), intent(in) :: STEPZ
524
REAL(dp), intent(in) :: origin(3)
525
REAL, intent(in) :: F(0:nmesh(1)-1,0:nmesh(2)-1,0:nmesh(3)-1)
526
REAL, intent(in) :: ZVALUE
527
REAL, intent(out) :: f2d(0:nmesh(1)-1,0:nmesh(2)-1)
528
529
530
INTEGER IC, IP, IPM, IX, K1, K2, K3
531
REAL(dp) zm, z
532
REAL f_cur, f_up, fv
533
534
if (nmesh(3) == 1) then
535
f2d(:,:) = f(:,:,0)
536
else
537
! Loop on mesh points
538
DO K1 = 0,NMESH(1)-1
539
DO K2 = 0,NMESH(2)-1
540
! z-direction is scanned from top to bottom
541
DO K3 = NMESH(3)-1,0,-1
542
543
! Calculate Z coordinate of this point:
544
Z = stepz * K3 + origin(3)
545
546
IF (Z .LT. ZVALUE) THEN
547
548
! Linear interpolation to find the value of F at ZVALUE
549
550
f_cur = f(k1,k2,k3)
551
f_up = f(k1,k2,k3+1)
552
ZM = stepz * (K3+1) + origin(3)
553
554
FV = f_up + (f_cur-f_up) * (ZVALUE - ZM) / (Z - ZM)
555
f2d(k1,k2) = FV
409
556
410
GOTO 10
411
412
ENDIF
413
414
ENDDO
415
416
WRITE(6,*) 'Z = ',ZVALUE,
417
. ' not found. It is probably outside your cell'
418
STOP
419
420
10 ENDDO
421
write(2,*)
422
ENDDO
423
424
CLOSE(2)
425
426
END
557
GOTO 10
558
559
ENDIF
560
561
ENDDO
562
563
WRITE(6,*) 'Z = ',ZVALUE, &
564
' not found. It is probably outside your cell'
565
STOP
566
10 ENDDO
567
ENDDO
568
569
endif
570
571
END SUBROUTINE ISOZ
572
Loggerhead is a web-based interface for Breezy
Version: 2.0.1