Editing Dos extract.py

=Instructions=  

* To use this script: 
 1) download it:  [[Image:dos_extract.tar.tgz]]
 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: (Please do not copy and paste the cotent to a .py file to run, Use the download link above)
 #!/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]]