Vtotav.py: Difference between revisions

From Wiki
Jump to navigation Jump to search
New page: #!/usr/bin/env python """ A script which averages a CHGCAR or LOCPOT file in one direction to make a 1D curve. User must specify filename and direction on command line. Depends on ase """ ...
 
No edit summary
 
(25 intermediate revisions by 6 users not shown)
Line 1: Line 1:
#!/usr/bin/env python
There is a script written in Python to convert LOCPOT for visualization.  
"""
A script which averages a CHGCAR or LOCPOT file in one direction to make a 1D curve.
User must specify filename and direction on command line.
Depends on ase
"""


import os
The source of this script is from Github [https://github.com/compphys/ase_tools/blob/master/scripts/vtotav.py] It is based on ASE; therefore you must first install the ASE environment either by using pip:
import sys
import numpy as np
import math
import string
import datetime
import time
from ase.calculators.vasp import VaspChargeDensity


starttime = time.clock()
  pip install --upgrade --user ase
print "Starting calculation at",
print time.strftime("%H:%M:%S on %a %d %b %Y")


if len(sys.argv) != 3:
or following the instructions on the webpage [https://wiki.fysik.dtu.dk/ase/install.html#installation-using-pip].
    print "\n** ERROR: Must specify name of file and direction on command line."
    print "eg. vtotav.py LOCPOT z."                                           
    sys.exit(0)


if not os.path.isfile(sys.argv[1]):
    print "\n** ERROR: Input file %s was not found." % sys.argv[1]
    sys.exit(0)


# Read information from command line
* 1) Download vtotav.py by clicking the link or from our Wiki: [[Image:Vtotav.tgz]]
# First specify location of LOCPOT
LOCPOTfile = sys.argv[1].lstrip()


# Next the direction to make average in
* 2) Uncompress vtotav.tgz and copy it to your bin folder ($HOME/bin)
# input should be x y z, or X Y Z. Default is Z.
allowed = "xyzXYZ"
direction = sys.argv[2].lstrip()
if allowed.find(direction) == -1 or len(direction)!=1 :
    print "** WARNING: The direction was input incorrectly." 
    print "** Setting to z-direction by default." 
if direction.islower():
    direction = direction.upper()
filesuffix = "_%s" % direction


# Open geometry and density class objects
* 3) Make it executable 
#-----------------------------------------
vasp_charge = VaspChargeDensity(filename = LOCPOTfile)
potl = vasp_charge.chg[-1]
atoms = vasp_charge.atoms[-1]
del vasp_charge


# For LOCPOT files we multiply by the volume to get back to eV
  chmod u+x ~/bin/vtotav.py
if 'LOCPOT' in LOCPOTfile:
    potl=potl*atoms.get_volume()


print "\nReading file: %s" % LOCPOTfile
* 4) In terminal, go to the folder containing work function calculation and run it with command:  
print "Performing average in %s direction" % direction


# Read in lattice parameters and scale factor
  vtotav.py  LOCPOT  z  (z is the direction, you can use x,y,and z, depending on your system)
#---------------------------------------------
cell = atoms.cell


# Find length of lattice vectors
* 5) You will get the output file named LOCPOT_Z. You can open it with Libreoffice/Excel or use one of the visualization methods given below.
#--------------------------------
latticelength = np.dot(cell, cell.T).diagonal()
latticelength = latticelength**0.5


# Read in potential data
#------------------------
ngridpts = np.array(potl.shape)
totgridpts = ngridpts.prod()
print "Potential stored on a %dx%dx%d grid" % (ngridpts[0],ngridpts[1],ngridpts[2])
print "Total number of points is %d" % totgridpts
print "Reading potential data from file...",
sys.stdout.flush()
print "done."


# Perform average
Vtotav.py is written in python 2. Depending on your version of the ASE environment it might not be compatible, as the newest ASE version requires python 3. Therefore you need to convert the python 2 script to python 3 by using the automated python 2 to 3 conversion programme [https://docs.python.org/2/library/2to3.html]:
#-----------------
if direction=="X":
    idir = 0
    a = 1
    b = 2
elif direction=="Y":
    a = 0
    idir = 1
    b = 2
else:
    a = 0
    b = 1
    idir = 2
a = (idir+1)%3
b = (idir+2)%3
# At each point, sum over other two indices
average = np.zeros(ngridpts[idir],np.float)
for ipt in range(ngridpts[idir]):
    if direction=="X":
        average[ipt] = potl[ipt,:,:].sum()
    elif direction=="Y":
        average[ipt] = potl[:,ipt,:].sum()
    else:
        average[ipt] = potl[:,:,ipt].sum()


if 'LOCPOT' in LOCPOTfile:
  2to3 -w vtotav.py
    # Scale by number of grid points in the plane.
    # The resulting unit will be eV.
    average /= ngridpts[a]*ngridpts[b]
else:
    # Scale by size of area element in the plane,
    # gives unit e/Ang. I.e. integrating the resulting
    # CHG_dir file should give the total charge.
    area = np.linalg.det([ (cell[a,a], cell[a,b] ),
                          (cell[b,a], cell[b,b])])
    dA = area/(ngridpts[a]*ngridpts[b])
    average *= dA


# Print out average
(-w to write it) and by changing the first line in the script to ''#!/usr/bin/env python3''
#-------------------
averagefile = LOCPOTfile + filesuffix
print "Writing averaged data to file %s..." % averagefile,
sys.stdout.flush()
outputfile = open(averagefile,"w")
if 'LOCPOT' in LOCPOTfile:
    outputfile.write("# Distance(Ang)    Potential(eV)\n")
else:
    outputfile.write("#  Distance(Ang)    Chg. density (e/Ang)\n")
xdiff = latticelength[idir]/float(ngridpts[idir]-1)
for i in range(ngridpts[idir]):
    x = i*xdiff
    outputfile.write("%15.8g %15.8g\n" % (x,average[i]))
outputfile.close()
print "done."


endtime = time.clock()  
 
runtime = endtime-starttime
 
print "\nEnd of calculation."
---------  Method 1 to visualize the LOCPOT_Z file
print "Program was running for %.2f seconds." % runtime
 
 
* Download this python script: [[Image:Workplot.tgz]] and in the terminal, run: python wplot.py
 
 
---------  Method 2 to visualize the LOCPOT_Z file
 
 
* You can also download another script to plot workfunction [[Image:Workplot_nathan.tgz‎]], but beware it gives huge figures :p (see below). The script requires one argument at the least; put "z" to define the file to be read as LOCPOT_Z
 
* It has additional options too; the first EXTRA argument passed in the commandline sets the title (use dashes ""). Add two additional numbers to specify the range on the x-axis.  
 
* Alternatively you can load the script into Excel, Libreoffice or Origin. You'll first have to convert the file to the right format by running the line below in the terminal:
 
  sed -i '1d' LOCPOT_Z &&  sed -i 's/  */\t/g' LOCPOT_Z &&  sed -i 's/^[ \t]*//' LOCPOT_Z
 
* [[Image: Workfunction.png]]

Latest revision as of 17:52, 10 December 2019

There is a script written in Python to convert LOCPOT for visualization.

The source of this script is from Github [1] It is based on ASE; therefore you must first install the ASE environment either by using pip: 

 pip install --upgrade --user ase

or following the instructions on the webpage [2].


  • 1) Download vtotav.py by clicking the link or from our Wiki: File:Vtotav.tgz
  • 2) Uncompress vtotav.tgz and copy it to your bin folder ($HOME/bin)
  • 3) Make it executable
 chmod u+x ~/bin/vtotav.py
  • 4) In terminal, go to the folder containing work function calculation and run it with command:
  vtotav.py  LOCPOT  z   (z is the direction, you can use x,y,and z, depending on your system)
  • 5) You will get the output file named LOCPOT_Z. You can open it with Libreoffice/Excel or use one of the visualization methods given below.


Vtotav.py is written in python 2. Depending on your version of the ASE environment it might not be compatible, as the newest ASE version requires python 3. Therefore you need to convert the python 2 script to python 3 by using the automated python 2 to 3 conversion programme [3]: 

  2to3 -w vtotav.py

(-w to write it) and by changing the first line in the script to #!/usr/bin/env python3


Method 1 to visualize the LOCPOT_Z file


  • Download this python script: File:Workplot.tgz and in the terminal, run: python wplot.py


Method 2 to visualize the LOCPOT_Z file


  • You can also download another script to plot workfunction File:Workplot nathan.tgz, but beware it gives huge figures :p (see below). The script requires one argument at the least; put "z" to define the file to be read as LOCPOT_Z
  • It has additional options too; the first EXTRA argument passed in the commandline sets the title (use dashes ""). Add two additional numbers to specify the range on the x-axis.
  • Alternatively you can load the script into Excel, Libreoffice or Origin. You'll first have to convert the file to the right format by running the line below in the terminal:
 sed -i '1d' LOCPOT_Z &&  sed -i 's/  */\t/g' LOCPOT_Z &&  sed -i 's/^[ \t]*//' LOCPOT_Z
Retrieved from "http://usrc.iciq.cat/wiki/index.php?title=Vtotav.py&oldid=7979"