-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdIdV_test_4N-coronene_cp2k.py
executable file
·198 lines (164 loc) · 9.39 KB
/
dIdV_test_4N-coronene_cp2k.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#!/usr/bin/python
import os
import numpy as np
import ppafm.io as io
import pyPPSTM as PS
import pyPPSTM.ReadSTM as RS
import matplotlib
# matplotlib.use('Agg') # Force matplotlib to not use any Xwindows backend. ## !!! important for working on clusters !!!!
import matplotlib.pyplot as plt
# --- specification of paths to the STM input files, PP positions and stored df results & format in which the data are stored - xsf or npy
path=''
path_pos='Q0.00K0.50/'
path_df = path_pos+'Amp0.40/'
data_format ="npy"
# --- specification of PBC, cell, and All the important stuff concerning electrons tunneling:
pbc=(0,0)
lvs = None
#lvs = np.array([[15., 0.,0.],[0.,15.,0],[0.,0.,15.]]
WorkFunction = 5.0 #more or less standart.
fermi=0.8933 # Normally is the Fermi Level read from the *.log file, but since CP2K gives the Fermi to the energy of HOMO and for this simulation we need it to be in the middle of the gap, we are changing it
orbs= 'sp' # 'sp' works now, 'spd' works for fireball as well
cut_min=-0.894 # HOMO 0.893 below the Fermi Level, other orbitals cut
cut_max=+0.894 # LUMO 0.893 eV above the Fermi Level
cut_at=-1 # All atoms of the molecule
eta = 0.01 # very low, to pronounce the single orbitals only
WF_decay=1.0 # for STM only - how fast the exponential decay fall, with the applied bias ( if 1 - 1:1 correspondence with bias; if 0, it doesn't change)
nV = 9 # for STM only - number of STM integrational steps nV ~ V/eta
lower_atoms=[] # No atoms has lowered hopping - be aware python numbering occurs here [0] - means lowering of the 1st atom
lower_coefs=[] # Lowering of the hoppings
# --- downloading and examples of downloading of the eigen-energies, the LCAO coefficients and geometry (this time for spin-unpolarized calculations):
#eigEn, coefs, Ratin = RS.read_FIREBALL_all(name = path+'phik_example_', geom=path+'crazy_mol.xyz', fermi=fermi, orbs = orbs, pbc=pbc,
# cut_min=cut_min, cut_max=cut_max,cut_at=cut_at, lower_atoms=lower_atoms, lower_coefs=lower_coefs);
#eigEn, coefs, Ratin = RS.read_AIMS_all(name = 'KS_eigenvectors_up.band_1.kpt_1.out', geom='geometry.in',fermi=fermi, orbs = 'sp', pbc=pbc,
# imaginary = False, cut_min=cut_min, cut_max=cut_max, cut_at=cut_at,
# lower_atoms=lower_atoms, lower_coefs=lower_coefs)
#eigEn, coefs, Ratin = RS.read_GPAW_all(name = 'out_LCAO_LDA.gpw', fermi=fermi, orbs = orbs, pbc=pbc,
# cut_min=cut_min, cut_max=cut_max, cut_at=cut_at, lower_atoms=lower_atoms, lower_coefs=lower_coefs);
eigEn, coefs, Ratin = RS.read_CP2K_all(name = 'crazy_mol', fermi=fermi, orbs = orbs, pbc=pbc,
cut_min=cut_min, cut_max=cut_max, cut_at=cut_at, lower_atoms=lower_atoms, lower_coefs=lower_coefs);
# --- the grid on which the STM signal is calculated; tip_r1 - PP distored by the relaxation in the PPAFM code; tip_r2 - uniform grid:
tip_r1, lvec, nDim, atomic_info_or_head = io.load_vec_field( path_pos+'PPpos' ,data_format=data_format)
dz=0.1
dx=dy =0.1
xl = lvec[1,0]
yl = lvec[2,1]
zl = lvec[3,2]
extent = (lvec[0,0],lvec[0,0]+xl,lvec[0,1],lvec[0,1]+yl)
tip_r2 = RS.mkSpaceGrid(lvec[0,0],lvec[0,0]+xl,dx,lvec[0,1],lvec[0,1]+yl,dy,lvec[0,2],lvec[0,2]+zl,dz)
# --- specification on which voltages the STM (dI/dV ...) calculations are performed - two methods - direct specification or sequence of voltages
Voltages=[-0.893,+0.894]
namez=['HOMO','LUMO']
#Voltages=np.arange(-1.0,+1.0+0.01,0.1) # this part is important for scans over slabs at different voltages
#namez = []
#for V in Voltages:
# namez.append(str(round(V,1)))
# --- downloading the df data
df, lvec2, nDim2, atomic_info_or_head = io.load_scal_field( path_df+'df' ,data_format=data_format)
# --- the Main Loop - for different WorkFunction (exponential z-decay of current), sample bias Voltages & eta - lorentzian FWHM
for WorkFunction in [WorkFunction]:
i=0;
for V in Voltages:
print ("Voltage:",V,"name:",namez[i])
for eta in [eta]:
current0 = PS.dIdV( V, WorkFunction, eta, eigEn, tip_r2, Ratin, coefs, orbs=orbs, s=1.0, px=0.0, py=0.0, pz = 0.0)
current1 = PS.dIdV( V, WorkFunction, eta, eigEn, tip_r1, Ratin, coefs, orbs=orbs, s=1.0, px=0.0, py=0.0, pz = 0.0)
current2 = PS.dIdV( V, WorkFunction, eta, eigEn, tip_r1, Ratin, coefs, orbs=orbs, s=0.0, px=1.0, py=1.0, pz = 0.0)
current3 = PS.dIdV_tilt( V, WorkFunction, eta, eigEn, tip_r1, tip_r2, Ratin, coefs, orbs=orbs, pz = 1.0, al =1.0)
current4 = PS.dIdV_tilt( V, WorkFunction, eta, eigEn, tip_r1, tip_r2, Ratin, coefs, orbs=orbs, dxyz = 1.0, al =1.0)
current5 = PS.STM( V, nV, WorkFunction, eta, eigEn, tip_r1, Ratin, coefs, orbs=orbs, px=0.5, py=0.5, WF_decay=WF_decay)
# next procedure is under development
current6 = PS.IETS_simple( V, WorkFunction, eta, eigEn, tip_r1, Ratin, coefs, orbs=orbs, s=0.0, px =0.5, py=0.5, pz=0.0, dxz=0.0, dyz=0.0, dz2=0.0, Amp=0.02)
# --- plotting part here, plots all calculated signals:
print(" plotting ")
for k in range(df.shape[0]):
dff = np.array(df[k,:,:]).copy()
curr0 = np.array(current0[k,:,:]).copy()
curr1 = np.array(current1[k,:,:]).copy()
curr2 = np.array(current2[k,:,:]).copy()
curr3 = np.array(current3[k,:,:]).copy()
curr4 = np.array(current4[k,:,:]).copy()
curr5 = np.array(current5[k,:,:]).copy()
curr6 = np.array(current6[k,:,:]).copy()
name_file='didV-'+namez[i]+'_%03d.dat' %k
name_plot_df='height:%03dA; df [Hz]' %k
name_plot0=namez[i]+';height:%03dA; dIdV [G0] s-fixed-tip' %k
name_plot1=namez[i]+';height:%03dA; dIdV [G0] s-tip' %k
name_plot2=namez[i]+';height:%03dA; dIdV [G0] pxy-tip' %k
name_plot3=namez[i]+';height:%03dA; dIdV [G0] pz-tip tilting' %k
name_plot4=namez[i]+';height:%03dA; dIdV [G0] dxyz-tip tilting' %k
name_plot5=namez[i]+';height:%03dA; STM [I] pxy-tip' %k
name_plot6=namez[i]+';height:%03dA; Under develoment' %k
# ploting part here:
plt.figure( figsize=(1.5* xl , 1.5*yl/2 ) )
plt.subplot(2,4,1)
plt.imshow( dff, origin='lower', extent=extent , cmap='gray')
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot_df)
# ploting part here:
plt.subplot(2,4,2)
plt.imshow( curr0, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot0)
# ploting part here:
plt.subplot(2,4,3)
plt.imshow( curr1, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot1)
# ploting part here:
plt.subplot(2,4,4)
plt.imshow( curr2, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot2)
plt.subplot(2,4,5)
plt.imshow( curr3, origin='lower', extent=extent , cmap='gray')
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot3)
# ploting part here:
plt.subplot(2,4,6)
plt.imshow( curr4, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot4)
# ploting part here:
plt.subplot(2,4,7)
plt.imshow( curr5, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot5)
# ploting part here:
plt.subplot(2,4,8)
plt.imshow( curr6, origin='lower', extent=extent, cmap='gray' )
plt.xlabel(r' Tip_x $\AA$')
plt.ylabel(r' Tip_y $\AA$')
plt.title(name_plot6)
plt.savefig( 'didv_cp2k_'+namez[i]+"_WF_"+str(WorkFunction)+"_"+str(eta)+'_%03d.png' %k , bbox_inches='tight' )
plt.close()
#plt.show()
#
# --- for saving WSxM format, if wanted
#
#tmp_curr=curr.flatten()
#out_curr=np.zeros((len(tmp_curr),3))
#out_curr[:,0]=tip_r[k,:,:,0].flatten()
#out_curr[:,1]=tip_r[k,:,:,1].flatten()
#out_curr[:,2]=tmp_curr.copy()
#f=open(name_file,'w')
#print >> f, "WSxM file copyright Nanotec Electronica"
#print >> f, "WSxM ASCII XYZ file; obtained from dIdV code by Krejci et al."
#print >> f, "X[A] Y[A] dIdV[G0]"
#print >> f, ""
#np.savetxt(f, out_curr)
#f.close()
#
#plt.show()
i = i+1
# --- the end
print()
print()
print("Done")