API tools faq
paste
Advertisement
Lucius_V

Partial DOS for ORCA

Lucius_V
Jun 3rd, 2015
265
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 8.48 KB | None | 0 0
raw download clone embed print report
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3.  
  4. """
  5. N-th molecular orbital represented as:
  6. \Psi(N) = \sum_m C_{mn}\phi_m
  7. where \phi_m is m-th atomic orbital
  8. It is sliced as Cmn[N,:]
  9. """
  10.  
  11. from pylab import *
  12. from scipy.stats import norm
  13. from scipy.constants import codata
  14. import sys
  15. import re
  16.  
  17. ATOMS_UNIQUE = False
  18. BROADENING = 0.1
  19. plevel={'NONE':0,'MINI':1,'NORMAL':2,'DEBUG':3}
  20. PRINTLEVEL = 'MINI'
  21. def gaussian(x, mu, sig):
  22.     return exp(-(x - mu)*(x - mu)/(2*sig*sig))
  23.  
  24. def list_to_array(List):
  25.     arr = []
  26.     for line in List:
  27.         for x in line[:-1].split():
  28.             arr.append(float(x))
  29.     return array(arr)
  30.  
  31. def plot_with_weight(value,weight):
  32.     gauss = norm(loc = value, scale = BROADENING)
  33.     domain = linspace(value-BROADENING*10,value+BROADENING*10)
  34.     plot(domain, gauss.pdf(domain))
  35.  
  36. def get_mo_decomposition(mo_number, data, renorm = False):
  37.     global DIM
  38.     global HFTyp
  39.     ln = 1
  40.     c = []
  41.     index = mo_number%6 + 2
  42.     shift = (DIM+4)*(mo_number/6)
  43.     for line in data:
  44.         if ln >= 5 + shift and ln <= (DIM + 4) + shift:
  45.             tmp = float(line[:-1].split()[index])
  46.             c.append(tmp)
  47.         ln += 1
  48.     arr = array(c)
  49.     if renorm and HFTyp == 'RHF':
  50.         return 2*arr*arr/(arr*arr).sum()
  51.     if renorm and HFTyp == 'UHF':
  52.         return arr*arr/(arr*arr).sum()
  53.     else:
  54.         return arr
  55.  
  56. def main():
  57.     global DIM
  58.     global HFTyp
  59.     for f in sys.argv[1:]:
  60.         if plevel[PRINTLEVEL] >= 1: print "Loading file "+f
  61.         species_found = False
  62.         species_amount_found = False
  63.         hftyp_found = False
  64.         nel_found = False
  65.         dim_found = False
  66.         mol_offset_found = False
  67.         with open(f,'r') as log:
  68.             ln = 1
  69.             data = []
  70.             data_a = []
  71.             data_b = []
  72.             species = []
  73.             for line in log:
  74.                 if not species_found:
  75.                     if "BASIS SET INFORMATION" in line:
  76.                         species_offset = ln
  77.                         species_found = True
  78.                 if species_found and ln == species_offset + 2:
  79.                     species_amount = int(line[:-1].split()[2])
  80.                     species_amount_found = True
  81.                 if species_amount_found and ln >= species_offset + 4 and ln <= species_offset + 4 + species_amount - 1:
  82.                     species.append(line[:-1].split()[3])
  83.                 if not hftyp_found:
  84.                     if "Hartree-Fock type" in line:
  85.                         HFTyp = str(line.split()[-1])
  86.                         if HFTyp == 'UKS':
  87.                             del(data)
  88.                         if HFTyp == 'RKS':
  89.                             del(data_a)
  90.                             del(data_b)
  91.                         hftyp_found = True
  92.                 if not nel_found:
  93.                     if "Number of Electrons" in line:
  94.                         NEL = int(line.split()[-1])
  95.                         nel_found = True
  96.                 if not dim_found:
  97.                     if "Basis Dimension" in line:
  98.                         DIM = int(line.split()[-1])
  99.                         dim_found = True
  100.                         if DIM%6 != 0: shift = 1
  101.                         else: shift = 0
  102.                 if not mol_offset_found:
  103.                     if "MOLECULAR ORBITALS" in line:
  104.                         MOL_OFFSET = ln
  105.                         REST = (MOL_OFFSET + 3)%212
  106.                         mol_offset_found = True
  107.                 if mol_offset_found and HFTyp == 'RHF' and ln >= MOL_OFFSET+2 and ln <= MOL_OFFSET+1+(DIM+4)*(DIM/6+shift):
  108.                     data.append(line[:-1])
  109.                 if mol_offset_found and HFTyp == 'UHF' and ln >= MOL_OFFSET+2 and ln <= MOL_OFFSET+1+(DIM+4)*(DIM/6+shift):
  110.                     data_a.append(line[:-1])
  111.                 if mol_offset_found and HFTyp == 'UHF' and ln >= MOL_OFFSET+2+(DIM+4)*(DIM/6+shift) +1 and ln <= MOL_OFFSET+1+2*(DIM+4)*(DIM/6+shift) +1:
  112.                     data_b.append(line[:-1])
  113.                 ln += 1
  114.             BASISAO = []
  115.             ATOMS_UNIQUE = False
  116.             if HFTyp == 'RHF':
  117.                 for line in data[4:4+DIM]:
  118.                     if not ATOMS_UNIQUE:
  119.                         BASISAO.append(re.sub("[0-9]","",line.split()[0]))
  120.                     else:
  121.                         BASISAO.append(line.split()[0])
  122.                         species = set(BASISAO)
  123.                         species = list(species)
  124.                         species.sort()
  125.                 if plevel[PRINTLEVEL] >= 2: print BASISAO
  126.                 if plevel[PRINTLEVEL] >= 1: print "Unique atoms:",species
  127.                 MOS_EIG = list_to_array(data[1::DIM+4])
  128.                 MOS_EIG = MOS_EIG*codata.value('Hartree energy in eV')
  129.                 MOS_OCC = list_to_array(data[2::DIM+4])
  130.                 domain = linspace(min(MOS_EIG)-1.0,max(MOS_EIG+1.0),10000)
  131.                 domain = linspace(-20,-4,10000)
  132.                 Cmn = []
  133.                 for n in arange(DIM):
  134.                     mo = get_mo_decomposition(n,data, True)
  135.                     Cmn.append(mo)
  136.                 Cmn = array(Cmn)
  137.                 BASISAO = array(BASISAO)
  138.                 for atom in species:
  139.                     Sum = zeros(domain.shape)
  140.                     indices = where(BASISAO == atom)
  141.                     for N,e in zip(arange(DIM),MOS_EIG):
  142.                         if plevel[PRINTLEVEL] >= 3: print e,Cmn[N,indices].sum()
  143.                         Sum += Cmn[N,indices].sum()*gaussian(domain,e,BROADENING)
  144.                     plot(domain,Sum,label=atom)
  145.                 Sum = zeros(domain.shape)
  146.                 for e in MOS_EIG:
  147.                     Sum += 2*gaussian(domain,e,BROADENING)
  148.                 plot(domain,Sum,'k',label='total DOS')
  149.             if HFTyp == 'UHF':
  150.                 for line in data_a[4:4+DIM]:
  151.                     if not ATOMS_UNIQUE:
  152.                         BASISAO.append(re.sub("[0-9]","",line.split()[0]))
  153.                     else:
  154.                         BASISAO.append(line.split()[0])
  155.                         species = set(BASISAO)
  156.                         species = list(species)
  157.                         species.sort()
  158.                 if plevel[PRINTLEVEL] >= 2: print BASISAO
  159.                 if plevel[PRINTLEVEL] >= 1: print "Unique atoms:",species
  160.                 MOS_EIG_A = list_to_array(data_a[1::DIM+4])
  161.                 MOS_EIG_A = MOS_EIG_A*codata.value('Hartree energy in eV')
  162.                 MOS_EIG_B = list_to_array(data_b[1::DIM+4])
  163.                 MOS_EIG_B = MOS_EIG_B*codata.value('Hartree energy in eV')
  164.                 MOS_OCC_A = list_to_array(data_a[2::DIM+4])
  165.                 MOS_OCC_B = list_to_array(data_b[2::DIM+4])
  166.                 domain = linspace(min(min(MOS_EIG_A),min(MOS_EIG_B))-1.0,max(max(MOS_EIG_A),max(MOS_EIG_B))+1.0,10000)
  167.                 domain = linspace(-20,-4,10000)
  168.                 Cmn_a = []
  169.                 Cmn_b = []
  170.                 for n in arange(DIM):
  171.                     mo_a = get_mo_decomposition(n,data_a, True)
  172.                     mo_b = get_mo_decomposition(n,data_b, True)
  173.                     Cmn_a.append(mo_a)
  174.                     Cmn_b.append(mo_b)
  175.                 Cmn_a = array(Cmn_a)
  176.                 Cmn_b = array(Cmn_b)
  177.                 BASISAO = array(BASISAO)
  178.                 for atom in species:
  179.                     Sum_a = zeros(domain.shape)
  180.                     Sum_b = zeros(domain.shape)
  181.                     indices = where(BASISAO == atom)
  182.                     for N,e in zip(arange(DIM),MOS_EIG_A):
  183.                         if plevel[PRINTLEVEL] >= 3: print e,Cmn_a[N,indices].sum()
  184.                         Sum_a += Cmn_a[N,indices].sum()*gaussian(domain,e,BROADENING)
  185.                     for N,e in zip(arange(DIM),MOS_EIG_B):
  186.                         if plevel[PRINTLEVEL] >= 3: print e,Cmn_b[N,indices].sum()
  187.                         Sum_b += Cmn_b[N,indices].sum()*gaussian(domain,e,BROADENING)
  188.                     plot(domain,Sum_a,label=atom)
  189.                     plot(domain,-Sum_b,label=atom)
  190.                 Sum_a = zeros(domain.shape)
  191.                 Sum_b = zeros(domain.shape)
  192.                 for e in MOS_EIG_A:
  193.                     Sum_a += gaussian(domain,e,BROADENING)
  194.                 for e in MOS_EIG_B:
  195.                     Sum_b += gaussian(domain,e,BROADENING)
  196.                 plot(domain,Sum_a,'k',label='total DOS')
  197.                 plot(domain,-Sum_b,'k',label='total DOS')
  198.             legend()
  199.             show()
  200.             log.close()
  201.  
  202. if __name__ == "__main__":
  203.     main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!