<jupyter_start><jupyter_text>Computing the nonadiabatic couplings in Kohn-Sham and excited states basesIn this tutorial, we will start computing the nonadiabatic couplings (NACs) from the molecular orbital overlap files obtained in [step2](../../7_step2_cp2k/1_DFT). The NACs will be computed in Kohn-Sham states and also in both single-partcile (SP) and many-body (MB) excited state bases. Finally we will plot the excited states energies vs time and NAC map. Table of contents1. [Importing needed libraries](import)2. [Overview of required files](required_files)3. [Computing the NACs](comp_nacs)         3.1. [Kohn-Sham basis](KS)\    3.2. [Excited state basis](excited_states)4. [Plotting the results](plotting)\    4.1. [NAC distribution](nac_dist)\    4.2. [Energy vs time](ene_time)\    4.3. [NAC map](nac_map)\    4.4. [Excitation analysis](exc_anal)\    4.5. [Average partial density of states](ave_pdos)      * 4.5.1. [Plot pDOS for all atoms angular momentums](ave_pdos_1)      * 4.5.2. [Plot pDOS for atoms with no angular momentum component](ave_pdos_2) A. Learning objectives* To be able to compute the NACs in Kohn-Sham and excited state bases* To be able to plot the NACs distribution * To be able to plot the computed excited states energies vs time* To be able to plot the NAC map* To be able to plot the average CI coefficients* To be able to plot the average partial density of states B. Use cases* [Computing the NACs](comp_nacs)* [Plotting the results](plotting) C. Functions- `libra_py`  - `CP2K_methods`    - [`read_cp2k_tddfpt_log_file`](exc_anal)  - `data_stat`    - [`cmat_distrib`](nac_dist)  - `workflows`    - `nbra`      - [`step2_many_body`](exc_anal)        - [`normalize_ci_coefficients`](exc_anal)      - [`step3`](comp_nacs)        - [`run_step3_ks_nacs_libint`](KS)        - [`run_step3_sd_nacs_libint`](excited_states)  - `units`    - `au2ev`  1. Importing needed libraries [Back to TOC](toc)Since the data are stored in sparse format using `scipy.sparse` library, we need to load this library so that we can read and check the orthonormality of the data.Import `numpy`, `scipy.sparse`, `CP2K_methods`, `data_stat`, `data_io`, `units`, `step2_many_body`, and `step3` modules. Also, `glob` will be needed to find specific types of files.<jupyter_code>import os
import sys
import time
import glob
import math
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from liblibra_core import *
from libra_py.workflows.nbra import step3, step2_many_body
from libra_py import CP2K_methods, units, data_stat, data_io<jupyter_output>/projects/academic/alexeyak/mohammad/software/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/alexeyak/mohammad/software/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/alexeyak/mohammad/software/miniconda3/envs/libra/lib/python3.6/import[...]<jupyter_text>2. Overview of required files [Back to TOC](toc)* `../../7_step2_cp2k/1_DFT/2_hpc/res`The MO overlap files are needed and stored in this folder.* `../../7_step2_cp2k/1_DFT/2_hpc/all_logfiles`All of the logfiles obtained from the electronic structure calculations of CP2K. These files will be needed to find the Kohn-Sham HOMO index and also extract the the TD-DFT data for each step.  3. Computing the NACs [Back to TOC](toc) 3.1. Kohn-Sham basis The `libra_py.workflow.nbra.step3.run_step3_ks_nacs_libint(params)` computes the NACs between pairs of Kohn-Sham states using the molecular orbital overlaps. The paramters for this function are as follows:`params['lowest_orbital']`: The lowest orbital considered in the computation of the MO overlaps. This value is exactly the sameas in the `run_template.py` file in step2.`params['highest_orbital']`: The highest orbital considered in the computation of the MO overlaps. This value is exactly the sameas in the `run_template.py` file in step2.`params['num_occ_states']`: The number of occupied orbitals to be considered from HOMO to lower occupied states. This value is defined by user.`params['num_unocc_states']`: The number of unoccupied orbitals to be considered from LUMO to higher unoccupied states. This value is defined by user.The two values above are used to create an active space which then will be used to select the elements from the MO overlap and energy matrices.`params['use_multiprocessing']`: A boolean flag to use the multiprocessing library of Python or not.`params['nprocs']`: The number of processors to be used for the calculations. Libra will use this only if the `params['use_multiprocessing']` is set to `True`.`params['time_step']`: The time-step used in the calculations in `fs`.`params['es_software']`: The name of the software package used to compute the electronic structure calculations. This will be used to generate the HOMO index of that system so it can build the active space.`params['path_to_npz_files']`: The full path to the MO overlap files.`params['logfile_directory']`: The full path to the folder where all the log files are stored.`params['path_to_save_ks_Hvibs']`: The full path to the folder in which the NACs between the Konh-Sham states are stored.`params['start_time']`: The start time-step.`params['finish_time']`: The finish time-step.`params['apply_phase_correction']`: A boolean flag for applying phase-correction algorithm.`params['apply_orthonormalization']`: A boolean flag for applying the orthonormalization algorithm.`params['do_state_reordering']`: If this value is set to `1` or `2`, the state-reordering will be applied to overlap matrices.`params['state_reordering_alpha']`: The state-reordering alpha value if the `params['do_state_reordering'] = 2`.After setting all the above paramters, the calculations are run using `step3.run_step3_ks_nacs_libint(params)`.<jupyter_code>params_ks = {
          'lowest_orbital': 24-10, 'highest_orbital': 24+11, 'num_occ_states': 10, 'num_unocc_states': 10,
          'use_multiprocessing': True, 'nprocs': 8, 'time_step': 1.0, 'es_software': 'cp2k',
          'path_to_npz_files': os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/res',
          'logfile_directory': os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/all_logfiles',
          'path_to_save_ks_Hvibs': os.getcwd()+'/res-ks-DFT',
          'start_time': 1200, 'finish_time': 1401,
          'apply_phase_correction': True, 'apply_orthonormalization': True,
          'do_state_reordering': 2, 'state_reordering_alpha':0
         }

#### For KS states - Applying correction to KS overlaps and computing the NACs in KS space
step3.run_step3_ks_nacs_libint(params_ks)<jupyter_output>Computing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in stepComputing orthonormalization of St matrices in step      1214 1221 1200120712281235

12491242





Done with step 1249 Done with step. Elapsed time:Done with step   1228Done with stepDone with stepDone with step0.20647907257080078Done with step 1242   
. Elapsed time: 1235Done with step 12001221Computing orthonormalization of St matrices in step1214 . Elapsed time:      1207 . Elapsed time:0.2748746871948242. Elapsed time:. Elapsed time:. Elapsed time:  12500.28510212898254395 

0.3064262866973877  0.3184823989868164. Elapsed time:


Computing orthonormalization of St matrices in step0.3124725818634033
 0.31462907791137695Computing ort[...]<jupyter_text>3.2. Excited state basis The NACs can also be computed between excited states. These include the single-particle and many-body bases which the latter is obtained from theTD-DFT calculations. First, we need to compute the overlap between excited state Slater-determinants (SDs) then they will be used to compute the NACsbetween them. For many-body states, the configuration interaction coefficietns will also be used. We will consider both single-particle and many-body for DFT calculations but only single-particle for xTB.To run the calculations `step3.run_step3_sd_nacs_libint(params)` function will be used. Some parameters are common with the ones used to run `step3.run_step3_ks_nacs_libint(params)`.There are different ways of defining the excited states SDs (the single-particle excited state basis). The first is through defining the `num_occ_states` and `num_unocc_states` in which Librawill start making the SDs from all of the occupied states (starting from `HOMO-num_occ_states+1`) to all of the unoccupied states (endsto `LUMO+num_unocc_states-1`). Also, if the unrestricted spin calculation flag is set to `True`, the SDs will be made for both alpha and beta spin channels. For example, if you want to build the electron-only excitation basis, you need to set `params['num_occ_states'] = 1` and set `params['num_unocc_states']`to a value less than the number of unoccupied orbitals that was considered in the computation of overlaps. This will generate all the electron-onlyexcitation from HOMO to unoccupied states.If the TD-DFT calculations has been done, then Libra will go over all log files and generate all the SDs used for all the steps and therefore the definition of these SDs is automatic and Libra will replace the `num_occ_states` and`num_unocc_states` itself based on the SDs that were generated from the TD-DFT log files. Other parameters needed to run the `step3.run_step3_sd_nacs_libint(params)` function are as follows:`params['isUKS']`: A boolean flag for unrestricted spin calculations.`params['is_many_body']`: If set to `True`, the NACs will be computed between pairs of many-body (TD-DFT) states. Also, the NACs between single-particle SDs obtained from the TD-DFT results will be computed as well. Otherwise, only single-particle NACs will be computed only for the SDs obtained from`num_occ_states` and `num_unocc_states`. This will be used for xTB calculations in which no TD-DFT is performed.`params['number_of_states']`: The number of TD-DFT states to consider. This value should not exceed the number of requested TD-DFT states in the CP2Kcalculations.`params['tolerance']`: A lower bound for selection of the excitation with configuration interaction coefficients higher than this value.`params['verbosity']`: An integer value showing the printing level. The default is set to 0. Higher values will print more data on terminal.`params['sorting_type']`: After defining the SDs, Libra will sort them either based on `'energy'` or `'identity'`.<jupyter_code>#### For excited states - Computing the excited states SDs and their overlaps and NACs
params_mb_sd = {
          'lowest_orbital': 24-10, 'highest_orbital': 24+11, 'num_occ_states': 10, 'num_unocc_states': 10,
          'isUKS': 0, 'number_of_states': 10, 'tolerance': 0.01, 'verbosity': 0,
          'use_multiprocessing': True, 'nprocs': 8,
          'is_many_body': True, 'time_step': 1.0, 'es_software': 'cp2k',
          'path_to_npz_files': os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/res',
          'logfile_directory': os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/all_logfiles',
          'path_to_save_sd_Hvibs': os.getcwd()+'/res-mb-sd-DFT',
          'outdir': os.getcwd()+'/res-mb-sd-DFT',
          'start_time': 1200, 'finish_time': 1401, 'sorting_type': 'identity',
          'apply_phase_correction': True, 'apply_orthonormalization': True,
          'do_state_reordering': 2, 'state_reordering_alpha':0
         }


step3.run_step3_sd_nacs_libint(params_mb_sd)<jupyter_output>sd_unique_basis is: [[[24, 25], 'alp'], [[24, 26], 'alp'], [[23, 25], 'alp'], [[22, 25], 'alp'], [[23, 26], 'alp'], [[22, 26], 'alp'], [[24, 27], 'alp'], [[24, 29], 'alp'], [[21, 26], 'alp'], [[20, 25], 'alp'], [[24, 28], 'alp'], [[20, 26], 'alp'], [[21, 25], 'alp'], [[21, 30], 'alp'], [[19, 26], 'alp'], [[19, 25], 'alp'], [[20, 27], 'alp'], [[23, 27], 'alp'], [[23, 28], 'alp'], [[24, 30], 'alp']]
sd_states_reindexed is: [[1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, -6], [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 7, -6], [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 8, -6], [1, -1, 2, -2, 3, -3, 4, -4, 7, -5, 6, -6], [1, -1, 2, -2, 3, -3, 7, -4, 5, -5, 6, -6], [1, -1, 2, -2, 3, -3, 4, -4, 8, -5, 6, -6], [1, -1, 2, -2, 3, -3, 8, -4, 5, -5, 6, -6], [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 9, -6], [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 11, -6], [1, -1, 2, -2, 8, -3, 4, -4, 5, -5, 6, -6], [1, -1, 7, -2, 3, -3, 4, -4, 5, -5, 6, -6], [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 10, -6], [1, -1, 8, -2, 3, -3, 4, -4, 5, -5, 6, -6], [...]<jupyter_text>4. Plotting the results [Back to TOC](toc) 4.1. NAC distribution One of the intuitive ways to visualize the NACs is to plot the distribution of the NACs. Here we plot them for SP and MB excited states.<jupyter_code>%matplotlib notebook

for basis in ['sd','ci']:
    nac = []
    nac_files = glob.glob(F'res-mb-sd-DFT/Hvib_{basis}*im*')
    for nac_file in nac_files:
        hvib = sp.load_npz(nac_file)
        hvib_dense = hvib.todense().real
        for i in range(hvib.shape[0]):
            for j in range(hvib.shape[0]):
                if j != i:
                    nac_ij = np.abs(hvib_dense[i,j])* 1000.0 * units.au2ev
                    x_mb = MATRIX(1,1)
                    x_mb.set(0, 0, nac_ij )
                    nac.append( x_mb )
    bin_supp, dens, cum = data_stat.cmat_distrib( nac, 0, 0, 0, 0, 50, 0.1)
    plt.plot( bin_supp, dens, label=basis)
plt.xlabel('|NAC|, meV')
plt.ylabel('PD, 1/meV')
plt.title('NAC distribution, SP vs MB')
plt.legend()
plt.tight_layout()
# plt.savefig('nac_dist_1.jpg', dpi=600)<jupyter_output><empty_output><jupyter_text>4.2. Energy vs time Here, we plot the excited states energy vs time. Since the excited states were sorted by their `'identity'` it is easy to visualize the states energies crossings.<jupyter_code>%matplotlib notebook
titles = ['SP basis', 'MB basis']
plt.figure()
for c, basis in enumerate(['sd','ci']):
    plt.subplot(1,2,c+1)
    energy_files = glob.glob(F'res-mb-sd-DFT/Hvib_{basis}*re*')
    energy_files = data_io.sort_hvib_file_names(energy_files)
    #print('Sorted energy files are:', energy_files)
    dt = 1.0 # fs
    energies = []
    for file in energy_files:
        energies.append(np.diag(sp.load_npz(file).todense().real))
    energies = np.array(energies)*units.au2ev
    md_time = np.arange(0,energies.shape[0]*dt,dt)
    for i in range(energies.shape[1]):
        plt.plot(md_time, energies[:,i]-energies[:,0])
    
    plt.title(titles[c])
    plt.ylabel('Excitation energy, eV')
    plt.xlabel('Time, fs')
    plt.tight_layout()<jupyter_output><empty_output><jupyter_text>4.3. NAC map Another way of visualizing the NAC values is to plot the average NAC matrix using `plt.imshow`.<jupyter_code>%matplotlib notebook
plt.figure()
nac_files = glob.glob('res-mb-sd-DFT/Hvib_sd*im*')
for c, nac_file in enumerate(nac_files):
    nac_mat = sp.load_npz(nac_file).todense().real
    if c==0:
        nac_ave = np.zeros(nac_mat.shape)
    nac_ave += np.abs(nac_mat)
nac_ave *= 1000*units.au2ev/c
nstates = nac_ave.shape[0]
plt.imshow(np.flipud(nac_ave), cmap='hot', extent=(0,nstates,0,nstates))#, vmin=0, vmax=100)
plt.xlabel('State index')
plt.ylabel('State index')
plt.colorbar().ax.set_title('meV')
plt.title('SP NACs')
%matplotlib notebook
nac_files = glob.glob('res-mb-sd-DFT/Hvib_ci*im*')
for c, nac_file in enumerate(nac_files):
    nac_mat = sp.load_npz(nac_file).todense().real
    if c==0:
        nac_ave = np.zeros(nac_mat.shape)
    nac_ave += np.abs(nac_mat)
nac_ave *= 1000*units.au2ev/c
nstates = nac_ave.shape[0]
plt.imshow(np.flipud(nac_ave), cmap='hot', extent=(0,nstates,0,nstates))#, vmin=0, vmax=100)
plt.xlabel('State index')
plt.ylabel('State index')
plt.colorbar().ax.set_title('meV')
plt.title('MB NACs')<jupyter_output><empty_output><jupyter_text>4.4. Excitation analysis Here, we average the CI coefficients for each TD-DFT excited states over the MD trajectory. This shows the excited states degree of mixing. This code is adopted from [here](https://github.com/AkimovLab/Project_CsPbI3_MB_vs_SP/blob/master/cubic/excitation_analysis/excitation_analysis.py) written by <NAME> with some modifications.We first need to define a dictionary which has exactly the same parameters as in `params_md_sd`. Other parameters that needs to be defined are as follows:`logfiles`: A list of all log files which is obtained from `glob.glob`.`nsds`: Number of SDs to be considered for each excited state.`colors`: The colors which are in order of the SDs. For example, the `green` color is for the first excitation that appears in the excited states analysis section in CP2K log files.The black lines show the error bars which are computed with a confidence interval coefficient of 1.96 which shows the confidence interval of $95\%$. The error bars are computed as follows:$$\epsilon=Z\frac{s}{\sqrt{N}}$$where $Z$ is the confidence interval coefficient, $s$ is the standard deviation of the $c_i^2$ values, and $N$ is the number of samples. We will use this formula over the tutorials for step4 as well.<jupyter_code>%matplotlib notebook
logfiles = glob.glob(F'{params_mb_sd["logfile_directory"]}/*.log')
params = params_mb_sd
nsds = 4
colors = ['green','blue','red','purple']
ci_coeffs = []
for logfile in logfiles:
    params.update({"logfile_name": logfile})
    excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = CP2K_methods.read_cp2k_tddfpt_log_file( params ) 
    ci_coefficients_raw_norm = step2_many_body.normalize_ci_coefficients(ci_coefficients_raw_unnorm)
    for j in range(len(ci_coefficients_raw_norm)):
        for k in range(len(ci_coefficients_raw_norm[j])):
            ci_coefficients_raw_norm[j][k] = ci_coefficients_raw_norm[j][k]**2
    ci_coeffs.append(ci_coefficients_raw_norm)


nsteps = len(ci_coeffs)
nstates = params["number_of_states"]

coeffs = []
coeffs_avg   = []
coeffs_error = []


for state in range(nstates):

    coeffs.append( [] )
    coeffs_avg.append( [] )
    coeffs_error.append( [] )

    for sd in range( nsds ):

        coeffs[state].append( [] )
        coeffs_avg[state].append( [] )
        coeffs_error[state].append( [] )

        for step in range( nsteps ):
            if len( ci_coeffs[step][state] ) < nsds and sd > len( ci_coeffs[step][state] )-1:
                coeffs[state][sd].append( 0.0 )
            else:
                coeffs[state][sd].append( ci_coeffs[step][state][sd] )
     
        mb_coeff_avg, mb_coeff_std = data_stat.scalar_stat( coeffs[state][sd] )
        coeffs_avg[state][sd].append( mb_coeff_avg )
        coeffs_error[state][sd].append( 1.96 * mb_coeff_std / math.sqrt(nsteps) )
  
        plt.plot(     state+1, mb_coeff_avg, color="black", marker='s', markerfacecolor=colors[sd], markeredgewidth=0.4, markersize=7)
        plt.errorbar( state+1, mb_coeff_avg, yerr=coeffs_error[state][sd], linestyle="None", color='black')


# plt.ylim(0,1)
plt.xlabel('State Index')
plt.ylabel('c$_{i}^2$')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>4.5. Average partial density of states In this section, we will plot the average partial density of states (pDOS) over the MD trajectory. There are two ways to take the average of the pDOS:1- Average all the pDOS files and then convolve the average pDOS for each element.2- Convolve the pDOS files and then take the average for each element.We choose the first one due to two reasons. First, the computational cost is much lower and we only need one convolution. Second is that averaging over the grid points (using the method 2) is dependent on the number of grid points we use for convolution which again adds to the complexity of the procedure. Here, we will use normalized Gaussian function for weighting the pDOS values and summing them.$$f(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$$This function is defined in the `gaussian_function` below. To apply this to a vector of numbers and sum all the weighted Gaussians, we use the `gaussian_function_vector` which will be used for pDOS plots.<jupyter_code>def gaussian_function(a, mu, sigma, num_points, x_min, x_max):
    pre_fact = (a/sigma)/(np.sqrt(2*np.pi))
    x = np.linspace(x_min, x_max, num_points)
    x_input = np.array((-1/2)/(np.square(sigma))*np.square(x-mu))
    gaussian_fun = pre_fact*np.exp(x_input)
    
    return x, gaussian_fun
    
def gaussian_function_vector(a_vec, mu_vec, sigma, num_points, x_min, x_max):
    for i in range(len(a_vec)):
        if i==0:
            sum_vec = np.zeros(num_points)
        energy_grid, conv_vec = gaussian_function(a_vec[i], mu_vec[i], sigma, num_points, x_min, x_max)
        sum_vec += conv_vec
    return energy_grid, sum_vec<jupyter_output><empty_output><jupyter_text>4.5.1. Plot pDOS for all atoms angular momentums In this part, we plot the pDOS for all of the angular momentum components of each atom. This is done by using the `orbital_cols`. In fact, the `orbital_cols` is related to `orbitals`. For example, for `s` orbital, we consider the 3rd index and for `p` orbital, we sum the columns from 4 to 6 (`range(4,7)`). Here we want to show how the code works and how the you can modify that based on your project. In the next section, we will show the pDOS only for atoms and sum all the components in each row of the pdos file. Other parameters are as follows:`atoms`: The atoms names which will be used in the labeling and plotting. The atoms order should be exactly the same as appear in the `.pdos` files. For example, the `*k1*.pdos` files contain the pDOS data for `Ti` atom and `*k2*.pdos` files contain the data for the `O` atom.`npoints`: The number of grid points for making the Gaussian functions. Note that, this value should be more than the number of states in the `.pdos` files.`sigma`: The standard deviation in eV.`shift`: This value shifts the minimum and maximum energy found in the `pdos_ave` and will extend the boundaries from both sides by `shift`eV.Finally, we will plot the total density of states. Note that the HOMO energy level is set to zero.<jupyter_code>%matplotlib notebook
path_to_all_pdos = os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/all_pdosfiles'
atoms = ['Ti', 'O']
orbitals_cols = [[3], range(4,7), range(7,12), range(12,19)]
orbitals = ['s','p','d','f']
npoints = 4000
sigma = 0.05 # eV
shift = 2.0 # eV
ave_pdos_convolved_all = []
for c1,i in enumerate([1,2]):
    pdos_files = glob.glob(path_to_all_pdos+F'/*k{i}*.pdos')
    for c2, pdos_file in enumerate(pdos_files):
        pdos_mat = np.loadtxt(pdos_file)
        if c2==0:
            pdos_ave = np.zeros(pdos_mat.shape)
        pdos_ave += pdos_mat
    pdos_ave /= c2+1
    pdos_ave[:,1] *= units.au2ev
    e_min = np.min(pdos_ave[:,1])-shift
    e_max = np.max(pdos_ave[:,1])+shift
    homo_level = np.max(np.where(pdos_ave[:,2]==2.0))
    homo_energy = pdos_ave[:,1][homo_level]
    for c3, orbital_cols in enumerate(orbitals_cols):
        try:
            sum_pdos_ave = np.sum(pdos_ave[:,orbital_cols],axis=1)
            ave_energy_grid, ave_pdos_convolved = gaussian_function_vector(sum_pdos_ave, pdos_ave[:,1], sigma,
                                                                               npoints, e_min, e_max)
            ave_pdos_convolved_all.append(ave_pdos_convolved)
            pdos_label = atoms[c1]+F', {orbitals[c3]}'
            plt.plot(ave_energy_grid-homo_energy, ave_pdos_convolved, label=pdos_label)
        except:
            pass


ave_pdos_convolved_total = np.sum(np.array(ave_pdos_convolved_all),axis=0)
plt.plot(ave_energy_grid-homo_energy, ave_pdos_convolved_total, color='black', label='Total')
plt.legend()
plt.xlim(-4,4)
plt.ylabel('pDOS, 1/eV')
plt.xlabel('Energy, eV')
plt.title('TiO$_2$, 300 K')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>4.5.2. Plot pDOS for atoms with no angular momentum component As you can see we have removed the `for` loop for the `orbital_cols` and in the `try` section we have set `sum_pdos_ave = np.sum(pdos_ave[:,3::],axis=1)` which will sum all the columns from 3rd index (`pdos_ave[:,3::]`).<jupyter_code>%matplotlib notebook
path_to_all_pdos = os.getcwd()+'/../../7_step2_cp2k/1_DFT/2_hpc/all_pdosfiles'
atoms = ['Ti', 'O']
npoints = 4000
sigma = 0.05
shift = 2.0 # eV
ave_pdos_convolved_all = []
for c1,i in enumerate([1,2]):
    pdos_files = glob.glob(path_to_all_pdos+F'/*k{i}*.pdos')
    for c2, pdos_file in enumerate(pdos_files):
        pdos_mat = np.loadtxt(pdos_file)
        if c2==0:
            pdos_ave = np.zeros(pdos_mat.shape)
        pdos_ave += pdos_mat
    pdos_ave /= c2+1
    pdos_ave[:,1] *= units.au2ev
    e_min = np.min(pdos_ave[:,1])-shift
    e_max = np.max(pdos_ave[:,1])+shift
    homo_level = np.max(np.where(pdos_ave[:,2]==2.0))
    homo_energy = pdos_ave[:,1][homo_level]
    try:
        sum_pdos_ave = np.sum(pdos_ave[:,3::],axis=1)
        ave_energy_grid, ave_pdos_convolved = gaussian_function_vector(sum_pdos_ave, pdos_ave[:,1], sigma,
                                                                           npoints, e_min, e_max)
        ave_pdos_convolved_all.append(ave_pdos_convolved)
        pdos_label = atoms[c1]
        plt.plot(ave_energy_grid-homo_energy, ave_pdos_convolved, label=pdos_label)
    except:
        pass


ave_pdos_convolved_total = np.sum(np.array(ave_pdos_convolved_all),axis=0)
plt.plot(ave_energy_grid-homo_energy, ave_pdos_convolved_total, color='black', label='Total')
plt.legend()
plt.xlim(-4,4)
plt.ylabel('pDOS, 1/eV')
plt.xlabel('Energy, eV')
plt.title('TiO$_2$, 300 K')
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12280235428137459, 'fraction_numerical': 0.03159293715587621, 'mean_word_length': 4.272472472472472, 'pattern_counts': {'":': 1, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '271081', 'license': 'CC0-1.0', 'n_tokens_mistral': 9667, 'n_tokens_neox': 8993, 'n_words': 2765, 'path': '6_dynamics/2_nbra_workflows/8_step3_cp2k/1_DFT/tutorials.ipynb', 'repo_name': 'langerest/Tutorials_Libra', 'size': '26342'}