Commit 3fc8a37a authored by Michal Hermanowicz's avatar Michal Hermanowicz
Browse files

Initial commit

parents
This diff is collapsed.
# Makefile for SURFACE
#
# Copyright (C) 2022 Michal Hermanowicz <mh@hermanowicz.eu>
#
# SURFACE is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# SURFACE is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with SURFACE. If not, see <https://www.gnu.org/licenses/>.
SHELL = /bin/sh
FC=gfortran
surface: surface.o getinput.o genraw.o analysis.o plot.o nms.o
$(FC) $^ -o $@
#libnms.a: libnms.o
# ar qc libnms.a $^
.PHONY: clean
clean:
$(RM) -rf surface *.o *.a
===================================================
SURFACE: Electronic surface states analysis program
===================================================
Copyright (C) 2022 Michal Hermanowicz <mh@hermanowicz.eu>
This file is part of SURFACE.
SURFACE is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
SURFACE is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with SURFACE. If not, see <https://www.gnu.org/licenses/>.
==============================
Third-party public domain code
==============================
The source file 'nms.f' contains third-party public domain code
written by several authors listed therein. The code, used for
interpolation and integration, is part of NMS package and comes from
the book:
D.K. Kahaner, C. Moler, and S. Nash. Numerical Methods and
Software. Prentice Hall, Englewood Cliffs, NJ, 1989.
See also: <https://gams.nist.gov/cgi-bin/serve.cgi/Package/NMS>
=================
About the program
=================
SURFACE is a post-processing program for electronic surface states
analysis based on wave functions extracted from DFT slab
calculations. It supports the Abinit software package
(https://www.abinit.org) and uses its 'cut3d' utility program to
analyse the output *_WFK files.
What it does:
- generates charge/probability density profiles as a function of the
lattice axis perpendicular to the crystal surface;
- interpolates and integrates the profiles within the desired bounds;
- applies a user-defined criterion to qualify the electronic state as
a surface or bulk state;
- outputs data ready for plotting with gnuplot.
Numerical wave functions, obtained from a DFT calculation, can be
analysed in the context of their spatial localisation which, for a
complex band structure, is usually a time consuming task. The Abinit's
*_WFK output file (written with 'prtwf' option set to 1) contains wave
function information at each band structure data point defined
uniquely by the k-point and energy pair. The Abinit's 'cut3d' utility
program enables one to extract real and imaginary parts of the wave
function along with their 3D indices within the calculated grid. With
some level of approximation, a surface state can be identified based
on spatial localisation of the charge (probability) density. The
criterion for such identification depends strongly on the system's
characteristics, e.g. slab thickness or chemical reactivity, and it
requires the user's consideration. See further below for papers
containing example results obtained that way.
=====================
Features and workflow
=====================
The program requires Abinit wave function file (*_WFK) generated in a
an Abinit slab calculation in order to perform surface/bulk states
analysis. Optionally (recommended) an *.agr band structure file can be
provided to obtain a color-coded band structure plot.
General workflow:
1. Extract real and imaginary parts of the wave functions from the
Abinit's *_WFK file,
2. Determine squared wave function module at each point,
3. Sum all the squared wave function modules on the grid's x-y planes
as a function of the z coordinate (z is perpendicular to the
surface). This results in the charge (probability) density profile
d(z),
4. Interpolate the d(z) profile function (with the functionality
provided in 'nms.f'),
5. Integrate the d(z) function with Cubic Hermite Function Integral
Evaluator ('nms.f'). This part depends on the user-provided
preconditions for the analysis. For example: If the system
contains 20 atomic monolayers, the user might define the
percentage of the total charge (probability) density that should
be localised within the first n monolayers in order to qualify the
given state as a surface state (or a bulk state otherwise).
6. Output the results, e.g.:
.wf_k106_b223_spinor1 BULK
.wf_k107_b223_spinor1 SS (1)
.wf_k108_b223_spinor1 SS (2)
The first column is a density profile filename containing the
k-point and band number that correspond to the band structure
calculated with Abinit; the second column is a resulting
identification which can be on of three: BULK, SS(1), SS(2), for
bulk and two possible surface states, respectively. The number in
brackets can be either 1 or 2 and it refers to the surface with
lower and higher z coordinate, respectively.
7. Optionally, output the gnuplot-formatted data ready for plotting
('results.gnuplot' file).
==========
Input file
==========
When started with no input file, the program will request information
from the standard input, e.g.:
SURFACE: Electronic surface states analysis program
Name of the Abinit *_WFK file: bi2se2o_WFK
nspinor=2 run? (1-Yes/0-No): 1
Spinor component (1-2): 1
Start at band number: 1
Stop at band number: 32
Start at K-point number: 1
Stop at K-point number: 216
Surface (1) integral limit A: 0.0
Surface (1) integral limit B: 29.53
Surface (2) integral limit A: 66.37
Surface (2) integral limit B: 94.0
OPTIONAL band structure *(.agr) file (0 if none): bi2se3o_EBANDS.agr
Silent mode? (Writing to files only) [y/n]: no
Most of the lines are rather self-explenatory, except the surface
integral limits -- those refer to the 'z' coordinates (perpendicular
to the surface) within the slab calculation supercell. Surface
integral limits should enclose several atomic monolayers of a given
slab surface, to which the identification criterion will be
applied. Currently, the criterion itself is hard-coded in 'analysis.f'
file (and can be easly changed): if 65% of the charge (probability)
density is localised on a given surface (between the limits A and B),
it is qualified as a surface state, and a bulk state otherwise.
The input data shown above can be provided as a file, e.g.:
surface < input.in | tee output.out
where 'input.in' contains:
bi2se2o_WFK
1
1
1
32
1
216
0.0
29.53
66.37
94.0
bi2se3o_EBANDS.agr
n
=====
Notes
=====
The program involves an intensive I/O. If a *.WFK file is very large,
it can be placed within a ramdisk to save the hardware I/O cycles.
The 'cut3d' utility ocassionaly makes an error in its output. It
should generate five columns of data: real and imaginary parts of the
wave function preceeded by the 3D indices. Sometimes a number appears
too long and two columns overlap making the line layout 4 columns
instead of 5. The SURFACE program handles it by skipping invalid files
and placing them in the 'errors.dat' file for later examination.
============
Requirements
============
The program is written in legacy Fortran. There might be an occasional
GNU extension used in the code, so GNU Fortran compiler is
recommended. Some versions of gfortran (especially newer ones) might
complain about the deprecated functions used or an obsolete fixed-form
source format -- both are harmless and can be neglected.
The Abinit's 'cut3d' utility program is required to extract data from
the *_WFK files.
============
Installation
============
To compile the program, simply type:
make
Note: It might be necessary to edit the Makefile to adjust it to your
environment. It has been tested with GNU Fortran 8.3.0.
It will create a single executable file called 'surface' in the
current directory which you can invoke as:
./surface
To make it accessible system-wide, simply move or copy it to the
directory included in your PATH variable.
=================
Scientific papers
=================
The method coded in SURFACE, in a different implementation, has been
used within published papers including the ones listed below. If you
found the program useful, you might want to consider citing one of
those as a reference.
[1] M. Hermanowicz, M.W. Radny; Iodine Adsorption on Bi2Se3
Topological Insulator; Physica Status Solidi - Rapid Research
Letters, 13(2), 1800460 (2019).
[2] M. Hermanowicz, M.W. Radny; Topological electronic states of
bismuth selenide thin films upon structural surface defects;
Comput. Mater. Sci. 117 (2016) 76-82.
\ No newline at end of file
C analysis.f - process electronic states extracted from *_WFK files
C
C Copyright (C) 2022 Michal Hermanowicz <mh@hermanowicz.eu>
C
C This file is part of SURFACE.
C
C SURFACE is free software: you can redistribute it and/or modify it
C under the terms of the GNU General Public License as published by
C the Free Software Foundation, either version 3 of the License, or
C (at your option) any later version.
C
C SURFACE is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with SURFACE. If not, see <https://www.gnu.org/licenses/>.
SUBROUTINE ANALYSIS (
* i, j, spinors, new, spinor_suffix, err,
* STATE, silent,
* band_start, k_start, band_stop, k_stop, datestr,
* s1_inta, s1_intb, s2_inta, s2_intb )
character(len=3) :: wf_file_kno
character(len=20) :: wf_file_bandno
character(len=25) :: wf_file
integer :: i, j, spinors, new, helper
integer :: readiostat
character(len=8) :: spinor_suffix
integer :: k_start, k_stop, band_start, band_stop, nlines,m,zets,z
double precision, allocatable, dimension(:) :: a,b,c,d,e,modsq
double precision, allocatable, dimension(:) :: mainResultz
double precision, allocatable, dimension(:) :: mainResultmod
double precision :: s, slabThick
double precision :: surf_1, surf_2
double precision :: s1_inta, s1_intb
double precision :: s2_inta, s2_intb
double precision, external :: dpchqa
double precision, allocatable, dimension(:) :: dd
integer :: err
double precision :: inta, intb, ss_depth1, ss_depth2
double precision :: full_integral, integral_s1, integral_s2
character(len=1) :: silent
CHARACTER*10 STATE
character(len=25) :: output_file
character(len=17) :: datestr
output_file = 'out_'//datestr//'.txt'
WRITE ( wf_file_kno, '(I0)') j
WRITE ( wf_file_bandno, '(I0)') i
wf_file = '.wf_k' // trim(adjustl(wf_file_kno)) // '_b' //
* trim(adjustl(wf_file_bandno)) // spinor_suffix
OPEN (unit = 10, file = trim ( wf_file ))
if ( i .eq. band_start .AND. j .eq. k_start ) then
OPEN (unit=17, file= output_file )
nlines = 0
do
read (10, *, END=2)
nlines = nlines + 1
end do
2 continue
end if
rewind 10
if ( .not. ALLOCATED(a) ) then
allocate (a(nlines))
end if
if ( .not. ALLOCATED(b) ) then
allocate (b(nlines))
end if
if ( .not. ALLOCATED(c) ) then
allocate (c(nlines))
end if
if ( .not. ALLOCATED(d) ) then
allocate (d(nlines))
end if
if ( .not. ALLOCATED(e) ) then
allocate (e(nlines))
end if
if ( .not. ALLOCATED(modsq) ) then
allocate (modsq(nlines))
end if
zets = 0
do m=1, nlines
read (10,*,iostat=readiostat) a(m), b(m), c(m), d(m),e(m)
if(readiostat /= 0) then
STATE = 'ERROR'
RETURN
end if
if ( c(m) .eq. c(m-1)) then
continue
else
zets = zets + 1
end if
end do
if (.not. ALLOCATED(mainResultz) ) then
allocate (mainResultz(zets))
end if
if (.not. ALLOCATED(mainResultmod) ) then
allocate (mainResultmod(zets))
end if
s=0
z=1
do m=1, nlines
modsq(m) = abs (cmplx ( d(m), e(m) ) )** 2
if ( c(m) .eq. c(m-1) .OR. m .eq. 1 ) then
s = s + modsq(m)
if (m .eq. nlines) then
C write(*,*) c(m), s
mainResultz(z) = c(m)
mainResultmod(z) = s
end if
else
C write(*,*) c(m-1), s
mainResultz(z) = c(m-1)
mainResultmod(z) = s
s=0
z=z+1
s=s+modsq(m)
end if
end do
if ( i .eq. band_start .AND. j .eq. k_start ) then
if ( s1_inta .lt. mainResultz(1) .or. s1_inta .gt.
* mainResultz(zets) ) then
write (*,*) ' Integral bound (s1_inta) out of range!'
stop
end if
if ( s1_intb .lt. mainResultz(1) .or. s1_intb .gt.
* mainResultz(zets) ) then
write (*,*) ' Integral bound (s1_intb) out of range!'
stop
end if
if ( s2_inta .lt. mainResultz(1) .or. s2_inta .gt.
* mainResultz(zets) ) then
write (*,*) ' Integral bound (s2_inta) out of range!'
stop
end if
if ( s2_intb .lt. mainResultz(1) .or. s2_intb .gt.
* mainResultz(zets) ) then
write (*,*) ' Integral bound (s2_intb) out of range!'
stop
end if
end if
if (.not. ALLOCATED(dd) ) then
allocate (dd(zets))
end if
call dpchez(size(mainResultz), mainResultz, mainResultmod,
* dd, .false., 0, 0, err)
full_integral = dpchqa(size(mainResultz), mainResultz,
* mainResultmod, dd, mainResultz(1), mainResultz(zets), err)
C write (*,*) 'full_integral = ', full_integral
integral_s1 = dpchqa(size(mainResultz), mainResultz,
* mainResultmod, dd, s1_inta, s1_intb, err)
C write (*,*) 'integral_s1', integral_s1
integral_s2 = dpchqa(size(mainResultz), mainResultz,
* mainResultmod, dd, s2_inta, s2_intb, err)
C write (*,*) 'integral_s2', integral_s2
C write (*,*)
if ( integral_s1 .ge. (0.65*full_integral) ) then
STATE = 'SS (1)'
else if ( integral_s2 .ge. (0.65*full_integral) ) then
STATE = 'SS (2)'
else
STATE = 'BULK'
end if
if ( silent .eq. 'n' ) then
write (17,*) wf_file, STATE
write (*,*) wf_file, STATE
else if ( silent .eq. 'y' ) then
write (17,*) wf_file, STATE
end if
CLOSE (10, status = 'DELETE')
if ( i .eq. band_stop .AND. j .eq. k_stop ) then
CLOSE (17, status = 'KEEP')
end if
END SUBROUTINE ANALYSIS
C genraw.f - generate raw data file with Abinit 'cut3d' tool
C
C Copyright (C) 2022 Michal Hermanowicz <mh@hermanowicz.eu>
C
C This file is part of SURFACE.
C
C SURFACE is free software: you can redistribute it and/or modify it
C under the terms of the GNU General Public License as published by
C the Free Software Foundation, either version 3 of the License, or
C (at your option) any later version.
C
C SURFACE is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with SURFACE. If not, see <https://www.gnu.org/licenses/>.
SUBROUTINE GENRAW ( WFK, i, j, nspinor2, spinors, spinor_suffix,
* band_start, k_start, band_stop, k_stop )
CHARACTER*59 run_cut3d / 'cut3d < .cut3d.in 1>/dev/null
* 2>./.cut3d.err' /
CHARACTER*20 WFK
INTEGER*4 status, system
CHARACTER*30 string / 'type cut3d >/dev/null 2>&1' /
CHARACTER*20 cut3d_command / 'which cut3d' /
character(len=8) :: spinor_suffix
integer :: i, j, nspinor2, spinors, new
integer :: k_start, band_start, band_stop, k_stop
if ( i .eq. band_start .AND. j .eq. k_start ) then
status = system ( string )
if ( status .ne. 0 ) then
write (*,*)
write (*,*) 'No cut3d tool found in your $PATH!'
write (*,*) ''
write (*,*) 'cut3d is distributed with the Abinit software'
write (*,*) 'package and is necessary for SURFACE to extract'
write (*,*) 'relevant data from the Abinit *_WFK files. You'
write (*,*) 'can download it from https://www.abinit.org'
write (*,*) ''
STOP
else
WRITE (*,"(A)", advance="no") ' Using cut3d binary: '
STATUS = SYSTEM ( cut3d_command )
end if
OPEN (unit = 9, file = ".cut3d.in", action="readwrite")
end if
write (unit = 9, fmt='(A)') WFK
write (unit = 9, fmt='(I3)') j
write (unit = 9, fmt='(I3)') i
if ( nspinor2 .eq. 1 ) then
write (unit = 9, fmt='(I1)') spinors
end if
write (unit = 9, fmt='(I1)') 0
write (unit = 9, fmt='(I1)') 0
write (unit = 9, fmt='(I1)') 4
write (unit = 9, fmt='(A)') '.wf'
write (unit = 9, fmt='(I1)') 0
C The following SYSTEM call starts cut3d to generate a single file.
status = system ( run_cut3d )
rewind 9
if ( i .eq. band_stop .AND. j .eq. k_stop ) then
CLOSE (9, status = 'KEEP')
end if
END SUBROUTINE GENRAW
C getinput.f - get user input
C
C Copyright (C) 2022 Michal Hermanowicz <mh@hermanowicz.eu>
C
C This file is part of SURFACE.
C
C SURFACE is free software: you can redistribute it and/or modify it
C under the terms of the GNU General Public License as published by
C the Free Software Foundation, either version 3 of the License, or
C (at your option) any later version.
C
C SURFACE is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License