<jupyter_start><jupyter_text>QE methods and QE_utilsIn this tutorial, we will explore various methods needed to handle Quantum Espresso (QE) calculations - to run them, prepare input, and extract output. All that will be done with the help of the **QE_methods** and **QE_utils** modules, which contains the following functions:**QE_methods*** cryst2cart(a1,a2,a3,r)* [Topic 2](topic-2) read_qe_schema(filename, verbose=0)* [Topic 3](topic-3) read_qe_index(filename, orb_list, verbose=0)* [Topic 4](topic-4) read_qe_wfc_info(filename, verbose=0)* [Topic 9](topic-9) read_qe_wfc_grid(filename, verbose=0)* [Topic 5](topic-5) read_qe_wfc(filename, orb_list, verbose=0)* read_md_data(filename)* read_md_data_xyz(filename, PT, dt)* read_md_data_xyz2(filename, PT)   * read_md_data_cell(filename)* out2inp(out_filename,templ_filename,wd,prefix,t0,tmax,dt)* out2pdb(out_filename,T,dt,pdb_prefix)* out2xyz(out_filename,T,dt,xyz_filename)* xyz2inp(out_filename,templ_filename,wd,prefix,t0,tmax,dt)* get_QE_normal_modes(filename, verbosity=0)* [Topic 1](topic-1) run_qe(params, t, dirname0, dirname1)                         * read_info(params)* read_all(params)* read_wfc_grid(params)**QE_utils*** get_value(params,key,default,typ)* split_orbitals_energies(C, E)* [Topic 7](topic-7) merge_orbitals(Ca, Cb)* post_process(coeff, ene, issoc)* [Topic 6](topic-6) orthogonalize_orbitals(C)* [Topic 8](topic-8) orthogonalize_orbitals2(Ca,Cb)<jupyter_code>import os
import sys
import math
import copy

if sys.platform=="cygwin":
    from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
    from liblibra_core import *
#from libra_py import *
from libra_py import units
from libra_py import QE_methods
from libra_py import QE_utils
from libra_py import scan

from libra_py import hpc_utils
from libra_py import data_read
from libra_py import data_outs
from libra_py import data_conv
from libra_py.workflows.nbra import step2


import py3Dmol   # molecular visualization
import matplotlib.pyplot as plt   # plots
%matplotlib inline 

plt.rc('axes', titlesize=24)      # fontsize of the axes title
plt.rc('axes', labelsize=20)      # fontsize of the x and y labels
plt.rc('legend', fontsize=20)     # legend fontsize
plt.rc('xtick', labelsize=16)    # fontsize of the tick labels
plt.rc('ytick', labelsize=16)    # fontsize of the tick labels

plt.rc('figure.subplot', left=0.2)
plt.rc('figure.subplot', right=0.95)
plt.rc('figure.subplot', bottom=0.13)
plt.rc('figure.subplot', top=0.88)

colors = {}

colors.update({"11": "#8b1a0e"})  # red       
colors.update({"12": "#FF4500"})  # orangered 
colors.update({"13": "#B22222"})  # firebrick 
colors.update({"14": "#DC143C"})  # crimson   

colors.update({"21": "#5e9c36"})  # green
colors.update({"22": "#006400"})  # darkgreen  
colors.update({"23": "#228B22"})  # forestgreen
colors.update({"24": "#808000"})  # olive      

colors.update({"31": "#8A2BE2"})  # blueviolet
colors.update({"32": "#00008B"})  # darkblue  

colors.update({"41": "#2F4F4F"})  # darkslategray

clrs_index = ["11", "21", "31", "41", "12", "22", "32", "13","23", "14", "24"]<jupyter_output>/projects/academic/cyberwksp21/Software/Conda/Miniconda3/envs/libra/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: to-Python converter for std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > already registered; second conversion method ignored.
  return f(*args, **kwds)
/projects/academic/cyberwksp21/Software/Conda/Miniconda3/envs/libra/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: to-Python converter for boost::python::detail::container_element<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, unsigned long, boost::python::detail::final_vector_derived_policies<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, false> > already registered; second conversion method ignored.
  return f(*args, **kwds)
/projects/academic/cyberwksp21/Software/Conda/Miniconda3/envs/libra/lib/python3.6/import[...]<jupyter_text>First, lets prepare the working directories and run simple SCF calculations to generate the output files<jupyter_code>PWSCF = os.environ['PWSCF62']
# Setup the calculations
params = {}

# I run the calculations on laptop, so no BATCH system
params["BATCH_SYSTEM"] = None

# The number of processors to use
params["NP"] = 1

# The QE executable
params["EXE"] = F"{PWSCF}/pw.x"

# The executable to generate the wavefunction files
params["EXE_EXPORT"] = F"{PWSCF}/pw_export.x" #"/mnt/c/cygwin/home/Alexey-user/Soft/espresso/bin/pw_export.x"

# The type of the calculations to be performed - in this case only a single SCF with spin-polarization
params["nac_method"] = 1
    
# The prefix of the input file
params["prefix0"] = "x0.scf" 

# Working directory - where all stuff happen
params["wd"] = os.getcwd()+"/wd"


# Remove the previous results and temporary working directory from the previous runs
os.system(F"rm -r {params['wd']}")
os.system(F"mkdir {params['wd']}")

# Copy the input files into the working directory
# also, notice how the SCF input file name has been changed
os.system(F"cp x0.scf.in {params['wd']}/x0.scf.0.in")
os.system(F"cp x0.exp.in {params['wd']}")
os.system(F"cp Li.pbe-sl-kjpaw_psl.1.0.0.UPF {params['wd']}")
os.system(F"cp H.pbe-rrkjus_psl.1.0.0.UPF {params['wd']}")<jupyter_output><empty_output><jupyter_text>1. run_qe(params, t, dirname0, dirname1)Use it to actually run the calculationsComment this out if you have already done the calculations<jupyter_code>help(QE_methods.run_qe)
!pwd
os.chdir("wd")
QE_methods.run_qe(params, 0, "res", "res2")
os.chdir("../")<jupyter_output>The time to run the QE calculations =  0.01<jupyter_text>2. read_qe_schema(filename, verbose=0)Can be used to read the information about the completed run<jupyter_code>pwd
info = QE_methods.read_qe_schema("wd/res/x0.save/data-file-schema.xml", verbose=0)
print(info)
nat = info["nat"]
R, F = info["coords"], info["forces"]

for at in range(nat):
    print(F"Atom {at} \t {info['atom_labels'][at]} \t\
          x={R.get(3*at+0):.5f}, y={R.get(3*at+1):.5f}, z={R.get(3*at+2):.5f}\
          fx={F.get(3*at+0):.5f}, fy={F.get(3*at+1):.5f}, fz={F.get(3*at+2):.5f}")<jupyter_output>Atom 0 	 Li 	          x=0.00000, y=0.00000, z=0.00000          fx=0.00000, fy=0.00000, fz=0.00000
Atom 1 	 H 	          x=2.83459, y=0.00000, z=0.00000          fx=0.00000, fy=0.00000, fz=0.00000<jupyter_text>3. read_qe_index(filename, orb_list, verbose=0)Is analogous to **read_qe_schema** in many regards, it just extracts a bit different info, including orbital energies.One would also need to specify which energy levels we want to extract, so one would need that info beforehands.In this example, we have just 4 electrons, so:1 - HOMO-12 - HOMO 3 - LUMO4 - LUMO+1Lets try just the 4 orbitals<jupyter_code>info2, all_e = QE_methods.read_qe_index("wd/res/x0.export/index.xml", [1,2,3,4], verbose=1)
print( info2)
print(all_e)
e_alp = all_e[0]
e_bet = all_e[1]
for i in range(4):
    print(F"E_{i}^alpha = {e_alp.get(i,i).real:12.8f} \t E_{i}^beta = {e_bet.get(i,i).real:12.8f}")<jupyter_output>E_0^alpha =  -1.86244053 	 E_0^beta =  -1.86219594
E_1^alpha =  -0.16462949 	 E_1^beta =  -0.16475681
E_2^alpha =  -0.08691250 	 E_2^beta =  -0.08650036
E_3^alpha =  -0.04025010 	 E_3^beta =  -0.04013592<jupyter_text>4. read_qe_wfc_info(filename, verbose=0)Can be used to extract some descriptors of the wavefunctions produced<jupyter_code>wfc_info1 = QE_methods.read_qe_wfc_info("wd/res/x0.export/wfc.1", verbose=1)
wfc_info2 = QE_methods.read_qe_wfc_info("wd/res/x0.export/wfc.2", verbose=1)
print(wfc_info1)
print(wfc_info2)<jupyter_output>{'ngw': 36889, 'igwx': 36889, 'nbnd': 10, 'nspin': 2, 'gamma_only': 'F', 'ik': 1, 'nk': 2}
{'ngw': 36889, 'igwx': 36889, 'nbnd': 10, 'nspin': 2, 'gamma_only': 'F', 'ik': 2, 'nk': 2}<jupyter_text>5. read_qe_wfc(filename, orb_list, verbose=0)Can be used to read in the actual wavefunctions produced<jupyter_code>alpha = QE_methods.read_qe_wfc("wd/res/x0.export/wfc.1", [1,2,3,4], verbose=0)
beta = QE_methods.read_qe_wfc("wd/res/x0.export/wfc.2", [1,2,3,4], verbose=0)
print(alpha)
print(alpha.num_of_rows, alpha.num_of_cols)

print(beta)
print(beta.num_of_rows, beta.num_of_cols)<jupyter_output><liblibra_core.CMATRIX object at 0x7fae6817b340>
36889 4
<liblibra_core.CMATRIX object at 0x7fae3ffb3730>
36889 4<jupyter_text>Orthogonality and normalizationBelow we can see that MO overlaps  are almost orthonormal - the diagonal elements are coorectly 1.0But the off-diagonal elements are not quite 0.0Same is true for  However, there is no any expectation about the orthogonality or normalization across the two sets<jupyter_code>S_aa = alpha.H() * alpha
S_bb = beta.H() * beta
S_ab = alpha.H() * beta
def print_mat(X):
    nr, nc = X.num_of_rows, X.num_of_cols
    for i in range(nr):
    
        line = ""
        for j in range(nc):
        
            line = line + "%8.5f  " % (X.get(i,j).real)
        print(line)

print("S_aa")
print_mat(S_aa)

print("S_bb")
print_mat(S_bb)

print("S_ab")
print_mat(S_ab)<jupyter_output>S_aa
 1.00000  -0.08150  -0.04761   0.03307  
-0.08150   1.00000   0.03712  -0.02467  
-0.04761   0.03712   1.00000  -0.01177  
 0.03307  -0.02467  -0.01177   1.00000  
S_bb
 1.00000  -0.08093  -0.04785   0.03336  
-0.08093   1.00000   0.03755  -0.02485  
-0.04785   0.03755   1.00000  -0.01246  
 0.03336  -0.02485  -0.01246   1.00000  
S_ab
 0.17301  -0.01049  -0.00412   0.00214  
-0.01635   0.15735   0.00452  -0.00236  
-0.00102  -0.00083  -0.06732   0.00139  
 0.01173  -0.00747  -0.00248   0.24420<jupyter_text>6. QE_utils.orthogonalize_orbitals(C)Can be used to orthogonalize orbitals if they are not.So lets transform alpha and beta orbitals such they are now orthonormal within each set.The resulting orbitals are not orthonormal across the two sets still<jupyter_code>alp = QE_utils.orthogonalize_orbitals(alpha)
bet = QE_utils.orthogonalize_orbitals(beta)

S_aa = alp.H() * alp
S_bb = bet.H() * bet
S_ab = alp.H() * bet

print("S_aa")
print_mat(S_aa)

print("S_bb")
print_mat(S_bb)

print("S_ab")
print_mat(S_ab)<jupyter_output>Testing if S is invertabile

[4, 1]
Det =  (0.9881907434095535-8.072456740942809e-20j)

Testing if S is invertabile

[4, 1]
Det =  (0.9883023952546621+8.762301957458487e-20j)
S_aa
 1.00000   0.00000   0.00000  -0.00000  
 0.00000   1.00000  -0.00000   0.00000  
 0.00000  -0.00000   1.00000   0.00000  
-0.00000   0.00000   0.00000   1.00000  
S_bb
 1.00000  -0.00000   0.00000   0.00000  
-0.00000   1.00000  -0.00000   0.00000  
 0.00000  -0.00000   1.00000  -0.00000  
 0.00000   0.00000  -0.00000   1.00000  
S_ab
 0.17301   0.00000  -0.00001   0.00000  
-0.00001   0.15735   0.00024  -0.00006  
-0.00001   0.00005  -0.06733   0.00028  
 0.00002   0.00017   0.00083   0.24420<jupyter_text>7. QE_utils.merge_orbitals(Ca, Cb)Sometimes (usually in the non-collinear case), we want to have a single set of orbitals (many are nearly doubly degenerate), not just alpha and beta components. We can prepare the single set from the spinor components using this function. In this example, we just gonna mimic non-collinear SOC calculations, pretending that alpha and beta orbital sets are the spinor components.<jupyter_code>C = QE_utils.merge_orbitals(alpha, beta)

S = C.H() * C
print_mat(S)<jupyter_output>1.00000  -0.08150  -0.04761   0.03307   0.17301  -0.01049  -0.00412   0.00214  
-0.08150   1.00000   0.03712  -0.02467  -0.01635   0.15735   0.00452  -0.00236  
-0.04761   0.03712   1.00000  -0.01177  -0.00102  -0.00083  -0.06732   0.00139  
 0.03307  -0.02467  -0.01177   1.00000   0.01173  -0.00747  -0.00248   0.24420  
 0.17301  -0.01635  -0.00102   0.01173   1.00000  -0.08093  -0.04785   0.03336  
-0.01049   0.15735  -0.00083  -0.00747  -0.08093   1.00000   0.03755  -0.02485  
-0.00412   0.00452  -0.06732  -0.00248  -0.04785   0.03755   1.00000  -0.01246  
 0.00214  -0.00236   0.00139   0.24420   0.03336  -0.02485  -0.01246   1.00000<jupyter_text>8. QE_utils.orthogonalize_orbitals2(Ca, Cb)This is a special orthogonalization procedure - the one for 2-component spinors. The inputs are assumed to be the components for each orbital. The orthogonalization works such that it is S_aa + S_bb = I<jupyter_code>alpha = QE_methods.read_qe_wfc("wd/res/x0.export/wfc.1", [1,2,3,4], verbose=0)
beta = QE_methods.read_qe_wfc("wd/res/x0.export/wfc.2", [1,2,3,4], verbose=0)
alp, bet = QE_utils.orthogonalize_orbitals2(alpha, beta)

S_aa = alp.H() * alp
S_bb = bet.H() * bet
print("S_aa")
print_mat(S_aa)

print("S_bb")
print_mat(S_bb)

print("S_aa + S_bb")
print_mat(S_aa + S_bb)

S_ab = alp.H() * bet
print("S_ab")
print_mat(S_ab)<jupyter_output>S_aa
 0.49999  -0.00014   0.00006  -0.00007  
-0.00014   0.49999  -0.00010   0.00004  
 0.00006  -0.00010   0.50000   0.00016  
-0.00007   0.00004   0.00016   0.50000  
S_bb
 0.50001   0.00014  -0.00006   0.00007  
 0.00014   0.50001   0.00010  -0.00004  
-0.00006   0.00010   0.50000  -0.00016  
 0.00007  -0.00004  -0.00016   0.50000  
S_aa + S_bb
 1.00000  -0.00000   0.00000  -0.00000  
-0.00000   1.00000  -0.00000   0.00000  
 0.00000  -0.00000   1.00000   0.00000  
-0.00000   0.00000   0.00000   1.00000  
S_ab
 0.08657  -0.00003   0.00001   0.00004  
-0.00005   0.07869   0.00008  -0.00006  
 0.00003  -0.00004  -0.03377   0.00018  
 0.00004   0.00006   0.00040   0.12212<jupyter_text>9. read_qe_wfc_grid(filename, verbose=0)Can be used to read the grid points for the given PW representation.<jupyter_code>G1 = QE_methods.read_qe_wfc_grid("wd/res/x0.export/grid.1", verbose=0)
print(len(G1))
for i in range(10):
    print(F"{i} \t {G1[i].x} \t {G1[i].y} \t {G1[i].z}")<jupyter_output>0 	 0.0 	 0.0 	 0.0
1 	 -1.0 	 0.0 	 0.0
2 	 1.0 	 0.0 	 0.0
3 	 -2.0 	 0.0 	 0.0
4 	 2.0 	 0.0 	 0.0
5 	 0.0 	 -1.0 	 0.0
6 	 0.0 	 0.0 	 -1.0
7 	 0.0 	 0.0 	 1.0
8 	 0.0 	 1.0 	 0.0
9 	 -1.0 	 -1.0 	 0.0
{'fraction_non_alphanumeric': 0.1457767614692011, 'fraction_numerical': 0.14384736315563812, 'mean_word_length': 3.926082365364308, 'pattern_counts': {'":': 13, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '3170', 'license': 'CC0-1.0', 'n_tokens_mistral': 7076, 'n_tokens_neox': 5851, 'n_words': 1278, 'path': '5_interactive_qe/1_LiH/test_QE_methods.ipynb', 'repo_name': 'compchem-cybertraining/Tutorials_QE_and_eQE', 'size': '13994'}