Removing duplicates PPSTM_simple

[ondrejkrejci]

Let's create a one function in pyPPSTM that would do anything from line 70 of PPSTM_simple.py;

• Adding atoms to the npy and xsf format for fixed scan to PPSTM_simple.py #32 can or cannot be part of the issue.

[LauriKurki]

This is a good idea. I was working on something along these lines a while back for the ASD-STM paper but never tested with other codes, fixed scans, output types etc. Heavily influenced by AFMulator in ppafm. The idea here is to create a STMulator object that can be initialised with scan parameters (params.ini, command line arguments),

PPSTM/pyPPSTM/STMulator.py

Lines 19 to 54 in c329d19

	class STMulator():
	"""
	Class for simulating scanning tunneling microscopy on molecules.
	Requires OpenCL_py3 branch of ProbeParticleModel to be installed to the same folder.
	Functionalities:
	1. Runs AFM scan on a molecule and saves the relaxed tip positions
	2. Runs relaxed STM/didv scan on a molecule with given parameters.
	3. Returns scan result as an np.ndarray
	Arguments:
	scan_type: string. STM or didv. (default: "didv")
	bias_voltage: float or tuple. (V or (V_min, V_max, dV))
	Bias voltage to use while scanning. Tuple format required for STM scan (defeult: 0.0)
	eta: float. broadening factor of the states. (default: 0.1)
	work_function: float. (eV) (default: 5.0)
	work_function_decay: float. rate of change for work function w.r.t. bias voltage (default: 0.5)
	tip_orb: string. the tip orbital [s, pxy, spxy, 5spxy, 10spxy, CO, pz, dz2, dxzyz] (default: s)
	sample_orb : string. orbitals of the sample 'sp' (light atoms only, faster) or 'spd' (all atoms) (default: 'sp')
	scan_window: tuple ((x1, y1, z1), (x2, y2, z2)). start and end positions for scanning in each direction.
	scan_dim: tuple (nx, ny, nz). number of tip steps in each direction.
	tip_type: string. [fixed, relaxed]. use relaxed or fixed tip. (default: 'relaxed')
	lvec: np.ndarray of shape (4, 3) or None. Unit cell boundaries for force field. First (row) vector
	specifies the origin, and the remaining three vectors specify the edge vectors of the unit cell.
	If None, will be calculated automatically from scan_window and tipR0, leaving some additional space
	on each side.
	pixPerAngstrome: int. Number of pixels (voxels) per angstrom in force field grid.
	iZPPs: int. Element of probe particle.
	qzs: Array of lenght 4. Position tip charges along z-axis relative to probe-particle center in angstroms.
	qs: Array of lenght 4. Values of tip charges in units of e/A^dim. Some can be zero.
	tip_r0: np.ndarray([x0, y0, z0]). Tip equilibrium position
	tip_stiffness: Array of lenght 4. [N/m] (x,y,z,R). stiffness of harmonic springs anchoring probe-particle to the metalic tip-apex
	relax_params: Array of lenght 4. (dt,damp, .., ..). parameters of relaxation, in the moment just dt and damp are used
	return_afm: boolean. Return afm scan of molecule also (default: False)
	timings: boolean. Print timing information. (default: False)
	"""

and then called with the molecule parameters

PPSTM/pyPPSTM/STMulator.py

Lines 121 to 135 in c329d19

	def eval(self, xyz, Zs, qs, eigs, coefs, REAs=None):
	"""
	Prepare and evaluate STM and AFM scan.
	Arguments:
	xyz: np.ndarray(N,3). Carteesian coordinates of each atom in the system
	Zs: np.ndarray(N). Atomic numbers of each atom
	qs: np.ndarray(N). Point charges of each atom
	eigs: np.ndarray. Eigenenergies of the system
	coefs: np.ndarray. KS coefficients of the system
	REAs: np.ndarray of shape (num_atoms, 4). Lennard Jones interaction parameters. Calculated automatically if None.
	Returns:
	if return_afm==True:
	tuple of np.ndarrays. STM/didv and AFM images
	else:
	np.ndarray. STM/didv image

Any thoughts on this approach?

[LauriKurki]

For this, I would also like to reorganise readSTM.py a bit, so that we pass cut_min, cut_max, cut_atoms etc. to the corresponding functions instead of using initial_check to set global variables. For example cut_eigenenergies would look something like

def cut_eigenenergies(eig, e_min, e_max):
  """
  Removes eigenenergies not within range [e_min, e_max].
  """

  mask = (eig > e_min) & (eig < e_max)
  return eig[mask]

[ondrejkrejci]

Ok sounds good to me, just I would do that in the beginning of the code so we spare memory . Sent from Outlook for Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Lauri Kurki ***@***.***> Sent: Tuesday, September 17, 2024 8:23:08 AM To: Probe-Particle/PPSTM ***@***.***> Cc: Krejci Ondrej ***@***.***>; Author ***@***.***> Subject: Re: [Probe-Particle/PPSTM] Removing duplicates PPSTM_simple (Issue #33) For this, I would also like to reorganise readSTM.py a bit, so that we pass cut_min, cut_max, cut_atoms etc. to the corresponding functions instead of using initial_check to set global variables. For example cut_eigenenergies would look something like def cut_eigenenergies(eig, e_min, e_max): """ Removes eigenenergies not within range [e_min, e_max]. """ mask = (eig > e_min) & (eig < e_max) return eig[mask] — Reply to this email directly, view it on GitHub<#33 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADFZAFLXGDXKXPBGEC6UVD3ZW64DZAVCNFSM6AAAAABN6PVLQ2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNJUGU2TAMZYGM>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[ondrejkrejci]

About the STMumlator - yes, that was exactly what I was thinking about. Is it working or not? If so, then the only thing is a missing documentation.

[LauriKurki]

It is working for FHI-aims input with relaxed scan and v-scan, but it needs to be tested for other settings before pushing to master branch.

[LauriKurki]

If we want to avoid the many files in tests, we could also probably switch to a configuration file parameter setup like in ppafm (using toml?).

My idea is that you could use

python run_ppstm.py config.toml

where config.toml contains the parameters we currently have on these lines in PPSTM_simple.py

PPSTM/PPSTM_simple.py

Lines 15 to 67 in 3e15f45

	ncpu = 1 # number of cpu cores for OMP paralelization: ncpu = 1 -- serial compilation & calculations; iff ncpu > 1, then OMP paralel recompilation is used and C++ calculations are running on more cores #
	#
	# ***** Main informations ******
	#
	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 = 'relaxed' # '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 #
	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 #
	WF_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 #
	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) #
	dft_code = 'fireball' # 'fireball'='Fireball'='FIREBALL' ; 'aims'='AIMS'='FHI-AIMS' ; 'cp2k'='CP2K' ; 'gpaw'='GPAW' #
	geometry_file = 'crazy_mol.xyz' # 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 = None # 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 = None # 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 #
	#
	# ***** Informations for x,y,z tip_type = 'fixed' ******
	#
	x = [ 0.0, 20.0, 0.25 ] # [xmin, xmax, dx] #
	y = [ 0.0, 15.0, 0.25 ] # [ymin, ymax, dy] #
	z = [ 10.0, 12.0, 0.1 ] # !!!! z-starts from zero - normally zmin >= 3-4 Ang above heighest atoms !!!! [zmin, zmax, dz] ; for single height scan use : [z, z, 1.0] #
	#
	# ***** Informations for PP positions, tip_type = 'relaxed' ******
	#
	Q = 0.00 # charge (PP-AFM) ; Ocharge PP-AFM complex_tip autumn 2018) ; [e] (monopole), [e*A] (dipole), [e*A^2] (quadrupole) #
	K = 0.50 # 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 #
	#
	# *****Output options ******
	#
	PNG = True # True / False -- plot "png" images (2D constant height) #
	WSxM = False # True / False -- write ".xyz" WSxM files (2D constant height) #
	XSF = 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 = False # 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 #
	plot_atoms = True # 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 #
	#
	# ***** Advanced options ******
	#
	cut_atoms = None # 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 ... #
	#
	# ***** More advanced options ******
	#
	WorkFunction = 5.0 # 5.0 eV is standart #
	fermi = None # None=0.0 -- no change to the Fermi Level ; -0.1 -- shifts the Fermi Level by 0.1 eV lower ... #
	cut_min = -2.5 # 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 = +2.5 # cut out all orbitals higher than -2.5 eV above Fermi (should be: cut_max >= Vmax+2*eta) . taken to the Fermi Level #
	files_path = '' # where are files fron DFT code ; rather do not use this #

With this approach it should be straightforward to add the cli scripts later for #34. What do you think?