Author Topic: Script for making data files to plot k-resolved transmission with gnuplot  (Read 28056 times)

0 Members and 1 Guest are viewing this topic.

Offline nori

  • QuantumATK Guru
  • ****
  • Posts: 122
  • Reputation: 12
    • View Profile
Dear all,

I'd like to share the script for calculating k-resolved transmission coefficients and making data files to plot them with gnuplot.

Converged 2probe NetCDF file is needed in order to use it.
In addition, you have to set three parameters in the script:
-line 16, NetCDF file name
-line 20, the number of k points
-line 23, the energy for which the transmission coefficients are to be calculated

Usage:
1. execute SCF calculation for your interesting system and save the result in NetCDF file
2. store "k_resolved_transmission.py" and the ncfile in the same directory
3. set-up three parameters mentioned above
4. execute "k_resolved_transmission.py"

After the task 4, you can get "xxxx.dat" for plotting.
(If the system is spin-polarized, "xxxx_up.dat" and "xxxx_dn.dat" are output.)

If you apply "k_resolved_transmission.py" to the system fe/mgo/fe (fe_mgo_fe.png),
you can obtain the picture like "fe_mgo_fe_parallel_0eV.GIF".
« Last Edit: February 5, 2009, 00:44 by nori »

Offline Nordland

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 812
  • Reputation: 18
    • View Profile
Very nice! I will try it out....

Offline ugglebot

  • Regular QuantumATK user
  • **
  • Posts: 22
  • Country: se
  • Reputation: 1
    • View Profile
Nori, could you please share also the GNUPlot plot commands you used to make the beautiful plot, from the data file?

Offline nori

  • QuantumATK Guru
  • ****
  • Posts: 122
  • Reputation: 12
    • View Profile
Hi ugglebot,

"pm3d" is used to plot "fe_mgo_fe_parallel_0eV.GIF".
So, If you get the error after executing command "set pm3d", please install the latest version of Gnuplot.

In the following, I show the commands:
gnuplot> set pm3d
gnuplot> set ticslevel 0
gnuplot> splot ‘xxxx.dat’ w pm3d

If you want 2D plot like "2d_plot.GIF", please input the commands as follows:
gnuplot> set pm3d map
gnuplot> set size square
gnuplot> splot ‘xxxx.dat’
« Last Edit: February 5, 2009, 01:09 by nori »

Offline Anders Blom

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 5541
  • Country: dk
  • Reputation: 91
    • View Profile
    • QuantumATK at Synopsys
Plotting k-resolved transmission with pylab directly from ATK
« Reply #4 on: February 13, 2009, 12:48 »
The ability to plot the transmission coefficients in various ways is an important addition to the basic functionality in VNL, so Nori has started a very relevant topic here! I would like to offer a second version, with focus on a few key points: The k-point grid. I hereby provide a Monkhorst-Pack class that can be used very generally. The main benefit is that the k-point generation is indepedent of the geometrical structure, so there is no need for complicated reading of the NetCDF file, and moreover there is a well-defined method to integrate over the k-points to obtain the total transmission. I will not document the Monkhorst-Pack class in this post; see the file, it should be quite self-explanatory. The plotting. I very much like Nori's plots, 3D plots are always nice. In case you are ok with the contour plot version only, it is however possible to do the plot with ATK itself, thus reducing the need for exporting the data to an external file, and use an external plotting program. Now, because the Monkhorst-Pack class is so general, the final user script becomes very small. Here is the part which does the calculations:
Code
from ATK.TwoProbe import *
from MonkhorstPack import MonkhorstPackGrid

nc_filename = 'femgofe_antipara.nc'
Npoints = 35
energy = 0.*eV

scf = restoreSelfConsistentCalculation(nc_filename)

kgrid = MonkhorstPackGrid((Npoints,Npoints),False)
K = kgrid.getFractionalKPoints()

trans_coeff = calculateTransmissionCoefficients(scf,energy,(K,Spin.Up))
Modify the code as needed for number of points, the energy, the NetCDF file name, and the spin component (the other one is Spin.Down). If you have an un-polarized calculation, replace the tuple (K,Spin.Up) with just plain K. Next, we prepare the data for plotting:
Code
Kxy = K[0:Npoints,1]
T = abs(trans_coeff.reshape(Npoints,Npoints))
Finally, the plotting part:
Code
import pylab as P
fig = P.figure()
P.xlabel('kx')
P.ylabel('ky')
P.contourf(Kxy,Kxy,T,40,cmap=P.cm.Spectral)
P.colorbar()
P.axis([-.5,.5,-.5,.5])
P.yticks(numpy.arange(-0.5,0.51,0.1))
P.xticks(numpy.arange(-0.5,0.51,0.1))
P.title('FeMgO, anti-parallel config, majority spin')
P.savefig('femgo_antipara_trans_up.png',dpi=100)
Play around with the plot parameters as you want (title, file name, etc)! The resulting plot of the example above is attached. (A side-comment: A subtle but hard-to-catch error will occur if you name the Monkhorst-Pack grid variable "grid", which seems like a natural thing to do, and then, instead of import pylab you think, "hey, why not do from pylab import * to make it simpler". The script will then fail, because the pylab import statement overwrites the name "grid" (it's a pylab function). This demonstrates the danger of from X import * statements; we generally recommend the use of import X instead (except I guess for from ATK.KohnSham import * :) !) The above script is simpler than Nori's (because I have hidden all the dirty stuff in the Monkhorst-Pack class), and therefore a bit easier to use. I think it is also a benefit to use a well-defined Monkhorst-Pack grid, for the reasons mentioned above. Note also that you don't really need the unit cell data to construct the k-point grid for the transmission, since we always use fractional coordinates for that. The disadvantage is that I did not consider time-reversal symmetry, as Nori did. So, mine is slower by a factor of 2, which probably is not be very critical in many cases, but it might be in some. I will, therefore, show in a separate post how we can take advantage of symmetries in the transmission spectrum, including also mirror planes etc, in a very powerful and time-consuming way. So, stay tuned - there is more to come!
« Last Edit: February 13, 2009, 12:51 by Anders Blom »

Offline M.Albert

  • Regular QuantumATK user
  • **
  • Posts: 13
  • Reputation: 0
    • View Profile
Dear nori
 I can obtain a very beautiful 3D picture with your method,but what physics properties can we analyze from this picture.
 By the way ,how to save the picture  in my computer .I am new to gnuplot ;D
THANK YOU!

Offline Nordland

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 812
  • Reputation: 18
    • View Profile
I am not a expert at analysis these pictures, but if you calculate the analytic result for the k-space resolved transmission for a tunneling barrier you will get the same picture as Nori shows in his post. Hence getting a k-space resolved transmission spectrum like this, will tell you that your scattering acts as a simple potential barrier.

Saving in gnuplot is a little tricky, but you can write something like:
set terminal png
set output "my_file.png"

And it should save the picture to my_file.png

Offline M.Albert

  • Regular QuantumATK user
  • **
  • Posts: 13
  • Reputation: 0
    • View Profile
Thank you ,I have tried it out.You are right!

Offline nori

  • QuantumATK Guru
  • ****
  • Posts: 122
  • Reputation: 12
    • View Profile
Dear M.Albert,

It might not be easy to take advantage of k-resolved transmission in order to investigate physical(transport) properties in the system,
but in some cases it's useful.

The system "fe/mgo/fe" is a perfect example.
If spin-configuration between electrodes is parallel, k-resolved transmission at fermi-level is like "fe_mgo_fe.png".
In "fe_mgo_fe.png", you can understand following things:
1. The nature of majority transmission is "Tunneling", because maximum value exists at gamma-point and transmission gets monotone decreasing when (kx, ky) gets larger.
(As Nordland mentioned above)
2. The nature of minority transmission is "Resonant Scattering" through interface state(s) between Fe and MgO, because there are some sharp peak away from gamma-point.

In addition, the nature of majority/minority transmission at fermi-level in the case of anti-parallel spin-configuration  is also "Resonant Scattering".
(Please see "femgo_antipara_trans_up.png" offerd by Anders)

These facts are very important for understanding TMR properties in the system "fe/mgo/fe".
For instance, if some adittional layers are inserted between Fe and MgO interface, It is expected that TMR ratio gets better.
Because interface states that contribute to "minority transmission in parallel" and "majority/minority transmission in anti-parallel" are broken by inserting addional layers, the conductances get dramatically decreasing while "majority transmission in parallel" gets little influence.
(Please refer to "Phys. Rev. B 72, 140404 (2005)" for more information)

Offline rosen

  • Heavy QuantumATK user
  • ***
  • Posts: 25
  • Reputation: 0
    • View Profile
Dear Anders Blom
If i want to plot the curve in an external program other than Pylab, how to pipe the output to a file such as XYZ?

Offline Anders Blom

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 5541
  • Country: dk
  • Reputation: 91
    • View Profile
    • QuantumATK at Synopsys
There is no standardized format for such files, so it depends a bit on which software you want to use. The data you need is kx and ky, which are contained in the variable Kxy (same grid in x and y), and the actual matrix T(kx,ky) which is the variable T. So, if your plotting software needs a format of the form
Quote
x1 y1 z x1 y2 z
you would loop over Kxy twice and print the matrix:
Code
for x in range(len(Kxy)):
    for y in range(len(Kxy)):
        print Kxy[x],Kxy[y],T[x,y]
    print
The extra print statement gives a new line after each sequence of y values. If this works or not depends , as mentioned, on the format expected. I guess you just have to play around with it, and try different things (obviously with a small number of k-points at first, not to waste time). Other software might just need the matrix itself and Kxy (I'm thinking of Matlab for instance). In that case print the matrix T and Kxy as they are, probably to separate files.

Offline kartiksamanta

  • New QuantumATK user
  • *
  • Posts: 3
  • Country: us
  • Reputation: 0
    • View Profile
Re: Plotting k-resolved transmission with pylab directly from ATK
« Reply #11 on: November 26, 2023, 02:56 »
Dear Anders Blom, I am using QuantumATK 2022.12 version. I was trying to use your script, plot_transmission.py, to extract the K-resolved 2D transmission but having the follwing error, "from MonkhorstPack import MonkhorstPackGrid ModuleNotFoundError: No module named 'MonkhorstPack' " Could you please modify the script for QuantumATK 2022.12 version to extract the k-resolved 2D transmission data and send to me! Thanking you! With best, Kartik Samanta ''''''''''''''''''''''''''''''''''''''''''''''''''''
The ability to plot the transmission coefficients in various ways is an important addition to the basic functionality in VNL, so Nori has started a very relevant topic here! I would like to offer a second version, with focus on a few key points: The k-point grid. I hereby provide a Monkhorst-Pack class that can be used very generally. The main benefit is that the k-point generation is indepedent of the geometrical structure, so there is no need for complicated reading of the NetCDF file, and moreover there is a well-defined method to integrate over the k-points to obtain the total transmission. I will not document the Monkhorst-Pack class in this post; see the file, it should be quite self-explanatory. The plotting. I very much like Nori's plots, 3D plots are always nice. In case you are ok with the contour plot version only, it is however possible to do the plot with ATK itself, thus reducing the need for exporting the data to an external file, and use an external plotting program. Now, because the Monkhorst-Pack class is so general, the final user script becomes very small. Here is the part which does the calculations:
Code
from ATK.TwoProbe import *
from MonkhorstPack import MonkhorstPackGrid

nc_filename = 'femgofe_antipara.nc'
Npoints = 35
energy = 0.*eV

scf = restoreSelfConsistentCalculation(nc_filename)

kgrid = MonkhorstPackGrid((Npoints,Npoints),False)
K = kgrid.getFractionalKPoints()

trans_coeff = calculateTransmissionCoefficients(scf,energy,(K,Spin.Up))
Modify the code as needed for number of points, the energy, the NetCDF file name, and the spin component (the other one is Spin.Down). If you have an un-polarized calculation, replace the tuple (K,Spin.Up) with just plain K. Next, we prepare the data for plotting:
Code
Kxy = K[0:Npoints,1]
T = abs(trans_coeff.reshape(Npoints,Npoints))
Finally, the plotting part:
Code
import pylab as P
fig = P.figure()
P.xlabel('kx')
P.ylabel('ky')
P.contourf(Kxy,Kxy,T,40,cmap=P.cm.Spectral)
P.colorbar()
P.axis([-.5,.5,-.5,.5])
P.yticks(numpy.arange(-0.5,0.51,0.1))
P.xticks(numpy.arange(-0.5,0.51,0.1))
P.title('FeMgO, anti-parallel config, majority spin')
P.savefig('femgo_antipara_trans_up.png',dpi=100)
Play around with the plot parameters as you want (title, file name, etc)! The resulting plot of the example above is attached. (A side-comment: A subtle but hard-to-catch error will occur if you name the Monkhorst-Pack grid variable "grid", which seems like a natural thing to do, and then, instead of import pylab you think, "hey, why not do from pylab import * to make it simpler". The script will then fail, because the pylab import statement overwrites the name "grid" (it's a pylab function). This demonstrates the danger of from X import * statements; we generally recommend the use of import X instead (except I guess for from ATK.KohnSham import * :) !) The above script is simpler than Nori's (because I have hidden all the dirty stuff in the Monkhorst-Pack class), and therefore a bit easier to use. I think it is also a benefit to use a well-defined Monkhorst-Pack grid, for the reasons mentioned above. Note also that you don't really need the unit cell data to construct the k-point grid for the transmission, since we always use fractional coordinates for that. The disadvantage is that I did not consider time-reversal symmetry, as Nori did. So, mine is slower by a factor of 2, which probably is not be very critical in many cases, but it might be in some. I will, therefore, show in a separate post how we can take advantage of symmetries in the transmission spectrum, including also mirror planes etc, in a very powerful and time-consuming way. So, stay tuned - there is more to come!

Offline Anders Blom

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 5541
  • Country: dk
  • Reputation: 91
    • View Profile
    • QuantumATK at Synopsys
This functionality is today included in the GUI, so I would recommend you use that instead. Anyway, it would not be possible to just port this script from such an old version, the API is completely different today.