init

Author: tzyii

ORCA_ipfitting.py

@@ -0,0 +1,227 @@
+#!/usr/bin/env python3
+
+def orca_ipfitting(orca_input, omega_l, omega_r, omega_tol=1.0e-3, maxiter=20):
+    import os, re, subprocess, shutil
+
+    inp_omega_decimal = 6
+    inp_omega_fmt = '%%.%df' % (inp_omega_decimal,)
+
+    omegas = []
+    types = []
+    E0s = []
+    E1s = []
+    IPs = []
+    kIPs = []
+
+    def inp_acc_omega(omega):
+        return round(omega, inp_omega_decimal)
+
+    def nearest_omega(omega, min_diff=None):
+        old_omega = None
+        for og in omegas:
+            diff = abs(og - omega)
+            if min_diff is None or diff < min_diff:
+                old_omega, min_diff = og, diff
+        return old_omega
+
+    def gen_fgbw_suffix(omega):
+        return '_' + inp_omega_fmt % (omega,)
+
+    def gen_inp(chg, spin, omega, finp):
+        with open(orca_input, 'r') as f:
+            txt = f.read()
+        xyz, xyzfile = True, True
+        txt_new, na = re.subn('\*\s*xyzfile\s+-?[0-9]+\s+-?[0-9]+', '*xyzfile %d %d' % (chg, spin), txt, re.I)
+        if na == 0:
+            txt_new, nb = re.subn('\*\s*xyz\s+-?[0-9]+\s+-?[0-9]+',     '*xyz %d %d' % (chg, spin), txt, re.I)
+            if nb == 0:
+                raise ValueError('Fail to renew charge and spin!')
+        txt_new, n = re.subn('RangeSepMu\s+[^\s]*', 'RangeSepMu   %.6f' % (omega), txt_new, re.I)
+        if n == 0:
+            raise ValueError('Fail to renew omega!')
+        with open(finp, 'w') as f:
+            f.write(txt_new)
+
+    def get_tote(flog):
+        tote = None
+        with open(flog, 'r') as f:
+            for line in f.readlines():
+                if line.startswith('FINAL SINGLE POINT ENERGY'):
+                    tote = float(line.strip().split()[-1])
+                    break
+        return tote
+
+    def get_homo(flog):
+        homo = None
+        with open(flog, 'r') as f:
+            lines = f.readlines()
+        start = None
+        for i, line in enumerate(lines):
+            if line.startswith('ORBITAL ENERGIES'):
+                start = i + 4
+                break
+        occ, lastocc, e, laste = 0.0, 0.0, 0.0, 0.0
+        for i, line in enumerate(lines[start:]):
+            occ, e = [float(x) for x in line.strip().split()[1:3]]
+            if round(occ) == 0 and round(lastocc) == 2:
+                homo = laste
+                break
+            lastocc, laste = occ, e
+        return homo
+
+    def calc_neutral(omega):
+        finp, flog, fgbw = 'neutral.inp', 'neutral.log', 'neutral.gbw'
+        gen_inp(0, 1, omega, finp)
+        old_omega = nearest_omega(omega)
+        if old_omega is not None:
+            old_fgbw = fgbw + gen_fgbw_suffix(old_omega)
+            shutil.copyfile(old_fgbw, fgbw)
+        with open(flog, 'w') as f:
+            exitcode = subprocess.call([orca_exec, finp], stdout=f, stderr=f)
+            if exitcode != 0:
+                raise RuntimeError('ORCA terminated abnormally!')
+        E, E_HOMO = get_tote(flog), get_homo(flog)
+        new_fgbw = fgbw + gen_fgbw_suffix(omega)
+        shutil.copyfile(fgbw, new_fgbw)
+
+        return E, E_HOMO
+
+    def calc_cation(omega):
+        finp, flog, fgbw = 'cation.inp', 'cation.log', 'cation.gbw'
+        gen_inp(1, 2, omega, finp)
+        old_omega = nearest_omega(omega, min_diff=0.05)
+        if old_omega is not None:
+            old_fgbw = fgbw + gen_fgbw_suffix(old_omega)
+        else:
+            old_fgbw = 'neutral.gbw'
+        shutil.copyfile(old_fgbw, fgbw)
+        with open(flog, 'w') as f:
+            exitcode = subprocess.call([orca_exec, finp], stdout=f, stderr=f)
+            if exitcode != 0:
+                raise RuntimeError('ORCA terminated abnormally!')
+
+        E = get_tote(flog)
+        new_fgbw = fgbw + gen_fgbw_suffix(omega)
+        shutil.copyfile(fgbw, new_fgbw)
+
+        return E
+
+    orca_root = os.environ.get('ORCA_ROOT', None)
+    if orca_root is None:
+        raise ValueError('"ORCA_ROOT" is not set!')
+    orca_exec = orca_root + '/orca'
+    if not os.access(orca_exec, os.X_OK):
+        raise ValueError('"%s" is not executable!')
+
+    E0r, E_HOMOr = calc_neutral(omega_r)
+    E1r = calc_cation(omega_r)
+
+    IPr, kIPr = E1r - E0r, -E_HOMOr
+    delta_r = IPr - kIPr
+
+    if IPr > kIPr:
+        raise ValueError("""\n***IP Fitting Error: Right Omega limit should have kIP > IP: {} !> {}""".format(kIPr, IPr))
+
+    omegas.append(omega_r)
+    types.append('Right Limit')
+    E0s.append(E0r)
+    E1s.append(E1r)
+    IPs.append(IPr)
+    kIPs.append(kIPr)
+
+    E0l, E_HOMOl = calc_neutral(omega_l)
+    E1l = calc_cation(omega_l)
+
+    IPl, kIPl = E1l - E0l, -E_HOMOl
+    delta_l = IPl - kIPl
+
+    if IPl < kIPl:
+        raise ValueError("""\n***IP Fitting Error: Left Omega limit should have kIP < IP: {} !< {}""".format(kIPl, IPl))
+
+    omegas.append(omega_l)
+    types.append('Left Limit')
+    E0s.append(E0l)
+    E1s.append(E1l)
+    IPs.append(IPl)
+    kIPs.append(kIPl)
+
+    converged = False
+    repeat_l = 0
+    repeat_r = 0
+    for step in range(maxiter):
+        # Regula Falsi (modified)
+        if repeat_l > 1:
+            delta_l /= 2.0
+        if repeat_r > 1:
+            delta_r /= 2.0
+        omega = - (omega_r - omega_l) / (delta_r - delta_l) * delta_l + omega_l
+        if inp_acc_omega(omega) == inp_acc_omega(omega_r) or inp_acc_omega(omega) == inp_acc_omega(omega_l):
+            converged = True
+            break
+
+        E0, E_HOMO = calc_neutral(omega)
+        E1 = calc_cation(omega)
+
+        IP, kIP = E1 - E0, -E_HOMO
+        delta = IP - kIP
+
+        if kIP > IP:
+            omega_r = omega
+            E0r = E0
+            E1r = E1
+            IPr = IP
+            kIPr = kIP
+            delta_r = delta
+            repeat_r = 0
+            repeat_l += 1
+        else:
+            omega_l = omega
+            E0l = E0
+            E1l = E1
+            IPl = IP
+            kIPl = kIP
+            delta_l = delta
+            repeat_l = 0
+            repeat_r += 1
+
+        omegas.append(omega)
+        types.append('Regula-Falsi')
+        E0s.append(E0)
+        E1s.append(E1)
+        IPs.append(IP)
+        kIPs.append(kIP)
+
+        # Termination
+        if abs(omega_l - omega_r) < omega_tol:
+            converged = True
+            break
+
+    print("""     => Regula Falsi Iterations <=\n""")
+    print("""    %3s %10s %14s %14s %14s %s""" % ('N','Omega','IP','kIP','Delta','Type'))
+    for k in range(len(omegas)):
+        print("""    %3d %10.6F %14.6E %14.6E %14.6E %s""" %
+               (k + 1, omegas[k], IPs[k], kIPs[k], IPs[k] - kIPs[k], types[k]))
+
+    if converged:
+        print("""\n    IP Fitting Converged""")
+        print("""    Final omega = %10.6F""" % ((omega_l + omega_r) / 2))
+
+    else:
+        raise ValueError("""IP Fitting """, step + 1)
+
+    return ((omega_l + omega_r) / 2)
+
+if __name__ == '__main__':
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--input', type=str, required=True, help='sample input file')
+    parser.add_argument('--min', type=float, required=True, help='minimum omega')
+    parser.add_argument('--max', type=float, required=True, help='maximum omega')
+    parser.add_argument('--tol', type=float, required=False, help='tolerance of omega', default=0.001)
+    parser.add_argument('--maxiter', type=int, required=False, help='max number of iterations', default=20)
+
+    args = parser.parse_args()
+
+    orca_ipfitting(args.input, args.min, args.max, args.tol, args.maxiter)
+

cry2cif.py

@@ -0,0 +1,73 @@
+#!/usr/bin/env python3
+
+from sys import argv
+import numpy as np
+
+def SkipLines(fd, nline):
+	for i in range(nline):
+		fd.readline()
+
+def WriteCif(param, atomsymbol, factcoor, filename = 'naphthalene'):
+	outp = []
+	atomamount = {}
+	atomformat = '   {0:<8s}1.0{1[0]:15.10f}{1[1]:15.10f}{1[2]:15.10f}   Biso{2:15.10f} {3:2s}\n'.format
+	outp.append('#' * 70 + '\n')
+	outp.append('#' + ' ' * 68 + '#\n')
+	outp.append('#' + ' ' * 23 + 'Optimized Crystal Data' + ' ' * 23 + '#\n')
+	outp.append('#' + ' ' * 68 + '#\n')
+	outp.append('#' + ' ' * 68 + '#\n')
+	outp.append('#' * 70 + '\n')
+	outp.append('\ndata_%s\n\n' % filename)
+	outp.append('%-38s%-15s\n' % ('_cell_length_a', param[0]))
+	outp.append('%-38s%-15s\n' % ('_cell_length_b', param[1]))
+	outp.append('%-38s%-15s\n' % ('_cell_length_c', param[2]))
+	outp.append('%-38s%-15s\n' % ('_cell_angle_alpha', param[3]))
+	outp.append('%-38s%-15s\n' % ('_cell_angle_beta', param[4]))
+	outp.append('%-38s%-15s\n' % ('_cell_angle_gamma', param[5]))
+	outp.append('%-38s%-15s\n' % ('_symmetry_space_group_name_H-M', "'P 1'"))
+	outp.append('%-38s%-15s\n\n' % ('_symmetry_Int_Tables_number', '1'))
+	outp.append("loop_\n_symmetry_equiv_pos_as_xyz\n   'x, y, z'\n\nloop_\n")
+	outp.append('   _atom_site_label\n   _atom_site_occupancy\n')
+	outp.append('   _atom_site_fract_x\n   _atom_site_fract_y\n   _atom_site_fract_z\n')
+	outp.append('   _atom_site_adp_type\n   _atom_site_B_iso_or_equiv\n   _atom_site_type_symbol\n')
+	for symbol, coor in zip(atomsymbol, factcoor):
+		atomamount[symbol] = atomamount.get(symbol, 0) + 1
+		atomname = symbol + str(atomamount[symbol])
+		outp.append(atomformat(atomname, coor, 1.0, symbol))
+	with open(filename + '.cif', 'w') as f:
+		f.write(''.join(outp))
+
+factcoor = []                  # Fractional coordinates
+lattice = []                   # Lattice vectors of supercell
+
+if len(argv) != 2:
+	print("Usage: cry2cif.py xx.out")
+elif argv[1][-4 :] == '.out':
+	readf = False
+	f = open(argv[1], 'r')
+	print('Reading ' + argv[1] + '...')
+	while True:
+		try:
+			line = f.readline()
+		except UnicodeDecodeError:
+			continue
+		if line[: 25] == ' FINAL OPTIMIZED GEOMETRY':
+			atomsymbol, factcoor = [], []
+			SkipLines(f, 5)
+			param = f.readline().split()       #lattice parameters
+			SkipLines(f, 4)
+			while True:
+				line = f.readline()
+				if line == '\n': break
+				sline = line.split()
+				atomsymbol.append(sline[3])
+				factcoor.append(sline[4 :])
+			readf = True
+			break
+	if not readf:
+		raise IOError('Final optimized geometry not found!')
+	f.close()
+	factcoor = np.array(factcoor, dtype = 'double')
+	WriteCif(param, atomsymbol, factcoor, argv[1].split('.')[0])     # Write cif for visualization
+else:
+    raise ValueError('Only output file of CRYSTAL17 named "xx.out" is supported!')

cryproperties.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+
+import os
+import re
+import numpy as np
+from io import StringIO
+from itertools import accumulate
+from functools import reduce
+from phonopy.units import Hartree
+
+for file in os.listdir(os.getcwd()):
+	if os.path.isfile(file):
+		if file[-4:] == '.out':
+			fout = file
+		elif file[-5:] == '.BAND':
+			fband = file
+		elif file[-5:] == '.DOSS':
+			fdoss = file
+
+with open(fout, 'r') as f:
+	outStr = f.read()
+kp = re.findall(r'LINE\s{0,}[0-9]+\s{0,}\((.*)\) IN TERMS OF PRIMITIVE LATTICE VECTORS\n\s{0,}([0-9]+)', outStr)
+paths = [re.sub(':', ' ---> ', p[0]) for p in kp]
+nums = [int(p[1]) for p in kp]
+accuNums = [0,] + list(accumulate(nums))
+efermi = max(float(e) for e in re.findall(r'TOP OF VALENCE BANDS.*EIG(.*)AU', outStr))
+
+with open(fband, 'r') as f:
+	bandStr = f.read()
+ticks = [float(t) for t in re.findall(r'@ XAXIS TICK\b.+,\s{0,}(.*)', bandStr)]
+lens = []
+for i in range(1, len(ticks)):
+	lens.append(ticks[i] - ticks[i - 1])
+band = np.loadtxt(StringIO(bandStr), comments=['#', '@'])[:, 1:] * Hartree
+
+print('%d paths as follow:' % len(paths))
+[print('    %s' % p) for p in paths]
+inp = input('please input the output order: ')
+order = [int(n) for n in inp.strip().split()]
+src = 0
+index, xaxis = [], []
+for n in order:
+	n_abs = abs(n)
+	index_tmp = list(range(accuNums[n_abs - 1], accuNums[n_abs]))
+	if n < 0:
+		index_tmp.reverse()
+	index += index_tmp
+	dst = src + lens[n_abs - 1]
+	xaxis.append(np.linspace(src, dst, nums[n_abs - 1]))
+	print('%9.6f ---> %9.6f' % (src, dst))
+	src = dst
+xaxis = np.column_stack([x.reshape(1, -1) for x in xaxis]).flatten()
+np.savetxt('band.dat', np.column_stack((xaxis, band[index, :])), fmt='%12.6f')
+
+doss = (np.loadtxt(fdoss, comments=['#', '@']) - np.array([efermi, 0])) * np.array([Hartree, 1 / Hartree])
+np.savetxt('doss.dat', doss, fmt='%12.6f')

daemon

@@ -0,0 +1,3 @@
+#!/usr/bin/env bash
+
+nohup ${*} &> /dev/null &

delpage.py

@@ -0,0 +1,17 @@
+#!/usr/bin/env python3
+
+from PyPDF2 import PdfFileReader, PdfFileWriter
+from sys import argv
+
+pageindex = int(argv[1])
+files = argv[2:]
+
+for file in files:
+	output = PdfFileWriter()
+	with open(file, 'rb') as infile:
+		input = PdfFileReader(infile)
+		for i in range(input.numPages):
+			if i != pageindex:
+				output.addPage(input.getPage(i))
+		with open('output-' + file, 'wb') as outfile:
+			output.write(outfile)