Dos extract.py
Instructions
- To use this script:
1) download it 2) copy it to the ~/bin folder 3) chmod u+x ~/bin/dos_extract.py 4) go to your calculation path which contains DOSCAR 5) run it with command in terminal: dos_extract.py
- The process to get the DOS informations
1) Read DOSCAR and get the NEDOS, atom numbers, etc in the calculations 2) Get the DOS information for specific atom 3) Get the DOS information for specific of the atom in 1) 4) Ask user to input the atoms 5) Ask user to input the orbitals 6) Sum the DOS for all the combinations.
- The content of this script:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
This script is used to get the specific DOS information.
For example: the 5th atoms' px orbital
Part of the codes are copied from the split_dos.py script of VTST.
A lots of ideas were also inspired by this script.
Author: Qiang Li
Email: lqcata@gmail.com or qli@bigbrosci.com
import numpy as np
from collections import defaultdict
### READ DOSCAR ###
def read_dosfile():
f = open("DOSCAR", 'r')
lines = f.readlines()
f.close()
index = 0
natoms = int(lines[index].strip().split()[0])
index = 5
nedos = int(lines[index].strip().split()[2])
efermi = float(lines[index].strip().split()[3])
return lines, natoms, nedos, efermi
### WRITE DOS0 CONTAINING TOTAL DOS ###
def write_dos0():
index = 5
lines, natoms, nedos, efermi = read_dosfile()
fdos = open("DOS0", 'w')
line = lines[index+1].strip().split()
ncols = int(len(line))
fdos.write('%15s,%15s,%15s \n' %('Energy', 'Spin-1', 'Spin-2'))
for n in xrange(nedos):
index +=1
e = float(lines[index].strip().split()[0])
e_f = e-efermi
spin_1 = float(lines[index].strip().split()[1])
spin_2 = float(lines[index].strip().split()[2])
fdos.write('%15.8f,%15.8f,%15.8f' % (e_f, spin_1, spin_2))
fdos.write('\n')
return index
### Get sprcific atoms dos information ###
def get_single_all(num):
index = 5
num = int(num)
atom_dos = []
lines, natoms, nedos, efermi = read_dosfile()
index = index + num * (nedos + 1 )
for n in xrange(nedos):
index +=1
#e = float(lines[index].strip().split()[0])
#e_f = e-efermi
dos = [float(i) for i in lines[index].strip().split()]
e = dos[0]
e_f = e -efermi
dos[0] = e_f
atom_dos.append(dos)
return atom_dos
orbital_nam = ['s','px','py','pz','dxy','dyz','dxz','dz2','dx2','f1','f2','f3','f4','f5','f6','f7']
orbital_num = range(1,17)
dict_orbital = {orbital_nam[i]:orbital_num[i] for i in range(0, len(orbital_nam))}
### Get sprcific orbital dos information for one specific atom ###
def get_single_orbital(num, orbital):
dos_all = get_single_all(num)
n_orbi = len(dos_all[0]) # n_orbi is the number of orbitals +1, this can be used to check ISPIN
atom_orbital = []
if n_orbi == 3 or n_orbi == 9 or n_orbi == 19 or n_orbi == 33:
for i in range(0, len(dos_all)):
dos_position = dict_orbital.get(orbital)
a, b, c = dos_all[i][0], dos_all[i][dos_position * 2 - 1], dos_all[i][dos_position * 2]
atom_orbital.append([a,b,c])
elif n_orbi == 2 or n_orbi == 5 or n_orbi == 10 or n_orbi == 17 :
for i in range(0, len(dos_all)):
dos_position = dict_orbital.get(orbital)
a, b = dos_all[i][0], dos_all[i][dos_position]
atom_orbital.append([a,b])
return atom_orbital
### READ CONTCAR file ###
def read_contcar():
f = open("CONTCAR", 'r')
lines = f.readlines()
f.close()
index = 5
ele_name = lines[index].strip().split()
index = 6
ele_num = [int(i) for i in lines[index].strip().split()]
dict_contcar = {ele_name[i]:ele_num[i] for i in range(0, len(ele_name))}
dict_contcar2 = defaultdict(list)
for element in ele_name:
indice = ele_name.index(element)
n_start = sum(ele_num[i] for i in range(0, indice+1)) - dict_contcar.get(element) +1
n_end = sum(ele_num[i] for i in range(0, indice+1)) +1
dict_contcar2[element].append(range(n_start, n_end))
return dict_contcar2
dict_contcar = read_contcar()
print '\nThe elements and their sequential numbers in CONTCAR are: \n'
for i in dict_contcar.keys():
list_for_usr = []
for j in dict_contcar.get(i)[0]:
list_for_usr.append(int(j))
print i, str(list_for_usr)
print '---' * 10 + '\n'
name = raw_input(Please Type the Atom(s) names or their sequential Numbers you are intersted in :
Choose 1) or 2) input style:
1) Use the element symbols:
Example A): Ru # For all atoms of sinlge element
Example B): Ru C H O # For all atoms of many elements
2) Use the sequential numbers in the POSCAR or CONTCAR:
Example A): 10 # for single atom: here is the 10th atom
Example B): 1 2 3 4 5 10 # for many atoms: here are the 1st, 2nd, 3rd, 4th, 5th and 10th atoms
Example C): 1-5 10 # the same results as B)
3) Note1): If the elements in your POSCAR or CONTCAR are like this: Ru C H O C (There are two C elements)
Please use the 2) input method. Otherwise only the first C atoms will be considered.
>>> )
ele_list = []
for i in name.split():
if '-' not in i:
if i.isdigit():
ele_list.append(int(i)) # in python the first one is 0
else : # get the numbers from element names
for j in read_contcar().get(i)[0]: # The first one is already 0.
ele_list.append(j)
else:
num_start = int(i.split('-')[0])
num_end = int(i.split('-')[1])+1
for j in range(num_start, num_end):
ele_list.append(j)
###--------------------------------------------###
orbital = raw_input( Please Type the Orbitals you are intersted in:
Choose 1) or 2) input style:
1) Use the orbital symbols:
Example A): p # for px, py and pz orbitals
Example B): d # for dxy, dyz, dxz, dz2 and dx2 orbitals
Example C): p d # sum of A) and B)
Example D): t2g # for dxy, dyz, dxz orbitals
Example E): eg # for dx2, dz2 orbitals
2) Use the specific orbital names:
Example A) px # to get the px orbital dos only
Example B) f1 # to get f1 orbital dos only
Example C) dxy dz2 # to get the sum dos of dxy and dz2 orbitals
Example D) px dxy f1 # to get the sum dos of px, dxy and f1 orbitals
3) Note: The d(x**2 - y**2) orbital is named dx2 here.
>>> )
### Creat a dictionary for 2) input style
label_1 = ['p', 'd', 'f', 't2g', 'eg']
orbital_1 =[ ['px','py','pz'], ['dxy','dyz','dxz','dz2','dx2'],
['f1','f2','f3','f4','f5','f6','f7'], ['dxy','dyz','dxz'], ['dz2','dx2']]
def convert_orbitals(orbital):
orbital = [dict_orbital.get(i) for i in orbital]
return orbital
dict_orbital2 = defaultdict(list)
for i in range(0, len(label_1)):
dict_orbital2[label_1[i]].append(orbital_1[i])
### Convert input information to orbital names ###
orbital_list = []
for i in orbital.split():
if i in label_1:
for j in dict_orbital2.get(i):
orbital_list.extend(j)
elif i in orbital_nam:
orbital_list.append(i)
### Convert the data into arrarys for sum
def convert_array(num, orbital):
test_data = get_single_orbital(num, orbital)
energy_l = []
orbital_dos_1_l = []
orbital_dos_2_l = []
if spin_or_not == 3 :
for i in range(0,len(test_data)):
energy_l.append(test_data[i][0])
orbital_dos_1_l.append(test_data[i][1])
orbital_dos_2_l.append(test_data[i][2])
energy = np.array(energy_l)
orbital_dos_1 = np.array(orbital_dos_1_l)
orbital_dos_2 = np.array(orbital_dos_2_l)
return energy, orbital_dos_1, orbital_dos_2
elif spin_or_not == 2:
for i in range(0,len(test_data)):
energy_l.append(test_data[i][0])
orbital_dos_1_l.append(test_data[i][1])
energy = np.array(energy_l)
orbital_dos_1 = np.array(orbital_dos_1_l)
return energy, orbital_dos_1
nedos = read_dosfile()[2]
spin_1 = np.zeros(nedos)
spin_2 = np.zeros(nedos)
### Chcek ISPIN ###
spin_or_not = len(get_single_orbital(2, 'px')[0])
if spin_or_not == 3:
print 'Spin unrestricted Calculation: ISPIN = 2 '
print 'Spin unrestricted Calculation: ISPIN = 2 '
elif spin_or_not == 2:
print 'Spin unrestricted Calculation: ISPIN = 1 '
print 'Spin unrestricted Calculation: ISPIN = 1 '
### Generate the sum dos data of orbitals for the atoms we selected ###
for i in ele_list:
for j in orbital_list:
if spin_or_not == 3 :
energy = convert_array(i, j)[0]
spin_1 = spin_1 + convert_array(i, j)[1]
spin_2 = spin_2 + convert_array(i, j)[2]
elif spin_or_not == 2:
energy = convert_array(i, j)[0]
spin_1 = spin_1 + convert_array(i, j)[1]
### Write the out-put file ###
out_name = raw_input(
Lets's give the out-put file a beautiful name:
For Example: Co-p.dat (Do not use Space in the name)
Type the out-put file name and click Enter:
>>> )
f_out = open(out_name.strip(), 'w')
for i in range(0, len(spin_1)):
if spin_or_not == 3 :
f_out.write('%15.8f %15.8f %15.8f \n' %(energy[i], spin_1[i], spin_2[i]))
elif spin_or_not == 2 :
f_out.write('%15.8f %15.8f \n' %(energy[i], spin_1[i]))
go back to Main Page, Computational Resources, Scripts, VASP, Núria López and Group