Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

33 ppstm simple to functions #36

Open
wants to merge 20 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 42 additions & 0 deletions configs/default.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
title = "Default configuration"

[scan]
scan_type = 'v-scan' # 'didv'='dIdV''='didv-single' -- only dIdV for one voltage = V ; 'v-scan'='V-scan'='Voltage-scan' -- both STM & dIdV scan - V .. Vmax; 'STM'='STM-single' -- STM for one Voltage = V, use V-scan rather; "states"="STATES" -- dIdV scancs at the eigenenergies, beware it still uses all other states for calculations #
tip_type = 'fixed' # 'fixed'='f' -- for stiff/metal tip apexes ; 'relaxed'='r' -- for flexible tip apexes (precalculated by PP-AFM) . For this option you have to have "installed" PPAFM in your PPSTM directory #
scan_window = [[0, 0, 10], [15.9, 15.9, 11.9]] # [[x0, y0, z0], [x1, y1, z1]] -- scan window in Angstroms, used for fixed scan #
scan_dim = [128, 128, 20] # [Nx, Ny, Nz] -- number of points in [x, y, z] used for fixed scan #
V = -0.5 # !!!! V = Vmin for SCAN !!!! #
V_max = +0.5 # V = V_min >= -2.0 V ; V_max <= 2.0 V (othervise changes in the later code needed) #
dV = 0.1 # # voltage step , dV <= 0.1 V #
eta = 0.1 # Lorentzian width of states in energy scale: typically 0.1; can be in range of 0.3-0.05 eV in some cases (low amount of layers ...) even up to 1.0 eV #
tip_orb = 's' # 's' ; 'pxy' -- px & py ; 'spxy' -- 50% s & 50% pxy ; '5spxy' -- 5% s & 95% pxy ; '10spxy' -- 10% s & 90% pxy ; 'CO' -- 13% s & 87% pxy (PRL 119, 166001 (2017)) ; 'pz' ; For sample_orbs = 'sp' , possible 'dz2' and 'dxzyz' -- dxz & dyz #
sample_orbs = 'sp' # orbitals of the sample 'sp' (light atoms only, faster) or 'spd' (all atoms) #

[input]
path = './' # where are files fron DFT code #
dft_code = 'aims' # 'fireball'='Fireball'='FIREBALL' ; 'aims'='AIMS'='FHI-AIMS' ; 'cp2k'='CP2K' ; 'gpaw'='GPAW' #
geometry_file = 'geometry.in' # E.G. 'input.xyz' , 'input.bas' , 'geometry.in'; None for GPAW #
pbc = [0,0] # (0,0) = None = False -- only original geometry ; (0.5,0.5) -- 2x2 cell ; (1,1) -- 3x3 cell (around original) ; (2,2) -- 5x5 cell (around original) ... #
lvs = [[1,0,0],[0,1,0],[0,0,1]] # None ; [[ax,ay,0],[bx,by,0]],[0,0,cz]] or [[ax,ay],[bx,by]] ; 'input.lvs' -- files with specified cell ; in FHI-AIMS & GPAW allready specified with geometry #
spin = false # None=False ; for FHI-AIMS & CP2K: None -- spin-unpolarized/spin-restricted calc. ; 'both' , 'up'='alpha' or 'down" (last 3 are for spin-polarizes or spin-unrestricted calculations) #
cp2k_name = 'crazy_mol' # Name used in CP2K calculations or GPAW calc #
Q = 0.0 # charge (PP-AFM) ; Ocharge PP-AFM complex_tip autumn 2018) ; [e] (monopole), [e*A] (dipole), [e*A^2] (quadrupole) #
K = 0.24 # x stiffness (PP-AFM master autumn 2018); klat (PP-AFM dev/OpenCl autumn 2018); Oklat (PP-AFM complex_tip autumn 2018) ; [N/m] #
data_format = 'npy' # 'xsf'='XSF' ; 'npy'='NPY' ; -- format in which PPpos are stored from PP-AFM run #

[advanced]
cut_atoms = -1 # None = -1 -- All atoms of the sample contributes to tunelling ; 1 -- only 1st atom of the sample contributes to the tunelling ; 57 -- first 57 atoms of the sample contributes to the tunelling ; ... #
lower_atoms = [] # [] = None -- No atoms has lowered hopping ; be aware python numbering occurs here: [0] - means lowering of the 1st atom; [0,1,2,3] -- lowering of 1st 4 atoms ... #
lower_coefs = [] # [] = None -- No lowering of the hoppings ; [0.5] -- lowering of the 1st atom hopping to 0.5 ; [0.5,0.5,0.5,0.5] -- lowering of 1st 4 atoms to 0.5 ... #
work_function = 5.0 # 5.0 eV is standart #
work_function_decay = 0.5 # 0.0 <= WF_decay <= 1.0 ; How fast WorkFunction tunnelling barrier is changing with Voltage : (WF = WF_0 + V*WF_decay) -- 0.0 no change ; 1.0 - the same change as voltage #
fermi = 0.0 # None=0.0 -- no change to the Fermi Level ; -0.1 -- shifts the Fermi Level by 0.1 eV lower ... #
cut_min = -3 # cut out all orbitals lower than -2.5 eV bellow Rermi (should be: cut_min <= Vmin-2*eta) . taken to the Fermi Level #
cut_max = 3 # cut out all orbitals higher than -2.5 eV above Fermi (should be: cut_max >= Vmax+2*eta) . taken to the Fermi Level #

[output]
plot_atoms = false # True / False -- plot "png" images (2D constant height) #
PNG = false # True / False -- write ".xyz" WSxM files (2D constant height) #
WSxM = false # True / False -- write ".xsf" files with 3D stucks of data . For this option you have to have "installed" PPAFM in your PPSTM directory #
NPY = true # True / False -- write ".npy" numpy binary files with 3D stucks of data . For this option you have to have "installed" PPAFM in your PPSTM directory #
XSF = false # True / False -- plot geometry (position of atoms into the PNG images and into the XSF files). You have to have your geometry, which you want to plot in input_plot.xyz. This doesn't change the name of the output files #
92 changes: 92 additions & 0 deletions ppstm_run.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
import os
import sys
import tomli

from pyPPSTM import basUtils as bU
from pyPPSTM import ReadSTM
from pyPPSTM import STMutils
from pyPPSTM import visualization

def main(config: dict):
# ppafm needed for relaxed tip scans and npy+xsf output
tip_type = config['scan']['tip_type']
npy_or_xsf_output = config['output']['NPY'] or config['output']['XSF']
if tip_type in ['relaxed', 'r'] or npy_or_xsf_output:
try:
import ppafm.io as io
except ImportError:
raise ImportError("ppafm is required for relaxed tip scans and NPY/XSF output.")

# Set tip coefficients
tip_orb = config['scan']['tip_orb']
tip_coeffs = STMutils.get_tip_coefficients(tip_orb)
print(f"Set tip coefficients for {tip_orb}: {tip_coeffs}")

# Read eigenenergies and coefficients
eigenenergies, coefs, atoms = ReadSTM.read_dft(config)
print(f"Read eigenenergies and coefficients.")

# Set up atom plotting if enabled
if config['output']['plot_atoms']:
try:
geom_plot, _, _ = bU.loadAtoms(
os.path.join(
config['input']['path'],
'input_plot.xyz'
)
)
except FileNotFoundError:
geom_plot = None
print("WARNING: Atom plotting disabled due to missing input_plot.xyz file.")

# Get tip positions
(
tip_r,
tip_r0,
lvec,
extent,
atomic_head_or_info
) = STMutils.get_tip_positions(config)
print(f"Tip positions read for a {config['scan']['tip_type']} scan.")

# Run STM scan
current, didv = STMutils.run_stm_scan(
config,
eigenenergies,
coefs,
tip_r,
atoms
)
print(f"STM scan complete for scan type {config['scan']['scan_type']}.")

# Get voltages and names
voltages, names = visualization.get_voltages_and_names(config, eigenenergies)

# Output results
if config['output']['PNG']:
visualization.plot_png(config, current, didv, voltages, names, lvec, extent, geom_plot)
print(f"PNG output complete.")
if config['output']['WSxM']:
visualization.plot_wsxm(config, current, didv, voltages, names, tip_r0)
print(f"WSXM output complete.")
if config['output']['XSF']:
visualization.save_xsf(config, current, didv, voltages, names, geom_plot, lvec)
print(f"XSF output complete.")
if config['output']['NPY']:
visualization.save_npy(config, current, didv, voltages, names, lvec, atomic_head_or_info)
print(f"NPY output complete.")

print(f"Output finished, exiting.")

if __name__=='__main__':
# Get config file from command line
config_file = sys.argv[1]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please add here a helper message, if there is no argv, or more than 1, thank you!


# Load config file
with open(config_file, 'rb') as f:
config = tomli.load(f)

print(f"Loaded config from {config_file}")
print(f"Config: {config}")

main(config)
174 changes: 159 additions & 15 deletions pyPPSTM/ReadSTM.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@

import time

from typing import Tuple

# this library has functions for reading STM coefficients and make a grid for non-relaxed 3D scan

# global variables:
Expand Down Expand Up @@ -201,18 +203,16 @@ def to_fermi(eig, fermi, orig_fermi=0.0):
eig -= fermi
return eig;

def cut_eigenenergies(eig):
def cut_eigenenergies(eig, cut_min, cut_max):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While, I am definitely for this one, I do not see, that the you would remove the global definitions at line 102 and 103. Would it be possible to remove them? If not, please add a comment so they should be removed in the future.

'''
Removes eigenstates (molecular orbitals) that are far from the energy important for scanning
'''
j = 1
global n_min_ , n_max_
for i in eig:
n_min_ = j if (i < cut_min_ ) else n_min_
n_max_ = j if (i < cut_max_ ) else n_max_
j += 1
assert (n_min_ < n_max_), "no orbitals left for dI/dV"
return eig[n_min_:n_max_];
'''
global n_min_, n_max_
n_min_ = sum(eig < cut_min)
n_max_ = sum(eig < cut_max)

mask = (eig > cut_min) & (eig < cut_max)
return eig[mask]

# procedure for handling the coefficients:

Expand Down Expand Up @@ -307,7 +307,7 @@ def read_AIMS_all(name = 'KS_eigenvectors.band_1.kpt_1.out', geom='geometry.in',
eig[i] = float(pre_eig[i])
del pre_eig, tmp;
eig = to_fermi(eig, fermi, orig_fermi=0.0)
eig = cut_eigenenergies(eig)
eig = cut_eigenenergies(eig, cut_min, cut_max)
print("eigenenergies read")
# ****** TO BE REMOVED ********
"""
Expand Down Expand Up @@ -399,7 +399,7 @@ def read_GPAW_all(name = 'OUTPUT.gpw', fermi = None, orbs = 'sp', pbc=(1,1), ima
eig = calc.get_eigenvalues(kpt=0, spin=0, broadcast=True)
at_num = slab.get_atomic_numbers()
eig = to_fermi(eig, fermi, orig_fermi=calc.get_fermi_level())
eig = cut_eigenenergies(eig)
eig = cut_eigenenergies(eig, cut_min, cut_max)
print("eigen-energies read")
# obtaining the LCAO coefficients (automatically removed unwanted states - molecular orbitals - and atoms)
coef = np.zeros((n_max_-n_min_,num_at_,Ynum_))
Expand Down Expand Up @@ -475,7 +475,7 @@ def read_FIREBALL_all(name = 'phi_' , geom='answer.bas', fermi=None, orbs = 'sp'
assert (len(eig)==n_bands), "number of bands wrongly specified"
eig = to_fermi(eig, fermi, orig_fermi=float(pre_eig[2]))
del pre_eig;
eig = cut_eigenenergies(eig)
eig = cut_eigenenergies(eig, cut_min, cut_max)
print("eigen-energies read")

print(" loading the LCAO coefficients")
Expand Down Expand Up @@ -540,7 +540,7 @@ def read_CP2K_all(name, fermi=None, orbs='sp', pbc=(1,1), imaginary = False, cut

# select relevant MOs
eig = to_fermi(eig, fermi, orig_fermi=fermi_energy)
eig = cut_eigenenergies(eig)
eig = cut_eigenenergies(eig, cut_min, cut_max)

# copy coefficients of relevant MOs
coef = np.zeros((n_max_-n_min_,num_at_,Ynum_))
Expand Down Expand Up @@ -686,4 +686,148 @@ def read_cp2k_MO_file(fn, spin):
# done
return labels, evals, occs, evecs, fermi_energy

############## END OF LIBRARY ##################################

def read_dft(
config: dict
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Reads DFT input using the configuration dictionary.

Args:
config (dict): Configuration dictionary.

Returns:
Tuple[np.ndarray, np.ndarray, np.ndarray]: Eigenvalues, coefficients, and atomic positions.
"""

input_params = config['input']
scan_params = config['scan']
advanced_params = config['advanced']

path = input_params["path"]
geometry_file = input_params["geometry_file"]
dft_code = input_params["dft_code"]
cp2k_name = input_params["cp2k_name"]
spin = input_params["spin"]

common_params = {
'fermi': advanced_params['fermi'],
'pbc': input_params['pbc'],
'orbs': scan_params['sample_orbs'],
'cut_min': advanced_params['cut_min'],
'cut_max': advanced_params['cut_max'],
'cut_at': advanced_params['cut_atoms'],
'lower_atoms': advanced_params['lower_atoms'],
'lower_coefs': advanced_params['lower_coefs']
}

if input_params['dft_code'].lower() in ('fireball', 'cp2k'):
lvs = input_params['lvs']
pbc = input_params['pbc']
if isinstance(lvs, (list, tuple, np.ndarray)):
cell = np.array(
[
[lvs[0][0], lvs[0][1], 0.0],
[lvs[1][0], lvs[1][1], 0.0],
[0.0, 0.0, 99.9]
]
) if (len(lvs) == 2) else lvs
elif isinstance(lvs, (str)):
cell = np.loadtxt(lvs)
elif pbc == [0, 0]:
cell = np.zeros((3,3))
else:
raise ValueError("pbc requested, but lvs not specified")

if ((dft_code == 'fireball') or(dft_code == 'Fireball') or (dft_code == 'FIREBALL')):
eigEn, coefs, Ratin = read_FIREBALL_all(
name=path+'phik_0001_',
geom=path+geometry_file,
lvs=cell,
**common_params
)

elif ((dft_code == 'gpaw') or(dft_code == 'GPAW')):
eigEn, coefs, Ratin = read_GPAW_all(
name=path+cp2k_name+'.gpw',
**common_params
)

elif ((dft_code == 'aims') or(dft_code == 'AIMS') or (dft_code == 'FHI-AIMS')):
if ((spin == None) or (spin == False)):
name = 'KS_eigenvectors.band_1.kpt_1.out'
elif spin in ['up', 'alpha']:
name = 'KS_eigenvectors_up.band_1.kpt_1.out'
elif spin in ['down', 'beta', 'dn']:
name = 'KS_eigenvectors_dn.band_1.kpt_1.out'
elif spin == 'both':
name_up = 'KS_eigenvectors_up.band_1.kpt_1.out'
name_dn = 'KS_eigenvectors_dn.band_1.kpt_1.out'
else :
raise ValueError(f"Unknown spin: {spin}")

if spin != 'both':
eigEn, coefs, Ratin = read_AIMS_all(
name=path+name,
geom=path+geometry_file,
**common_params
)
else:
eigEn1, coefs1, Ratin = read_AIMS_all(
name=path+name_up,
geom=path+geometry_file,
**common_params
)
eigEn2, coefs2, Ratin = read_AIMS_all(
name=path+name_dn,
geom=path+geometry_file,
**common_params
)
eigEn = np.concatenate((eigEn1, eigEn2), axis=0)
coefs = np.concatenate((coefs1, coefs2), axis=0)

elif ((dft_code == 'cp2k') or(dft_code == 'CP2K')):
if ((spin == None)or(spin == False)):
eigEn, coefs, Ratin = read_CP2K_all(
name = path + cp2k_name ,
lvs=cell,
**common_params
)
elif ((spin == 'up')or(spin == 'alpha')):
eigEn, coefs, Ratin = read_CP2K_all(
name=path+cp2k_name,
lvs=cell,
spin='alpha',
**common_params
)
elif (spin == 'both'):
eigEn1, coefs1, Ratin = read_CP2K_all(
name=path+cp2k_name,
lvs=cell,
spin='alpha',
**common_params
)
eigEn2, coefs2, Ratin = read_CP2K_all(
name = path + cp2k_name ,
lvs=cell,
spin='beta',
**common_params
)
eigEn = np.concatenate((eigEn1, eigEn2), axis=0)
coefs = np.concatenate((coefs1, coefs2), axis=0)
elif ((spin == 'down')or(spin == 'beta')or(spin == 'dn')):
eigEn, coefs, Ratin = read_CP2K_all(
name=path+cp2k_name,
lvs=cell,
spin='beta',
**common_params
)
else :
raise ValueError(f"Unknown spin: {spin}")

else:
raise ValueError(f"Unknown DFT code: {dft_code}")

return eigEn, coefs, Ratin

############## END OF LIBRARY ##################################
Loading