#***********************************************************
# * Copyright (C) 2017-2019 Brendan A. Smith, Wei Li, and Alexey V. Akimov
# * This file is distributed under the terms of the
# * GNU General Public License as published by the
# * Free Software Foundation; either version 3 of the
# * License, or (at your option) any later version.
# * http://www.gnu.org/copyleft/gpl.txt
#***********************************************************/
"""
.. module:: step2_analysis
:platform: Unix, Windows
:synopsis: This module implements functions for postprocessing the data generated by the
step2 module functions
.. moduleauthor:: Brendan A. Smith, Wei Li, and Alexey V. Akimov
"""
import os
import sys
import math
from libra_py.workflows.nbra import step2_many_body
# Fisrt, we add the location of the library to test to the PYTHON path
if sys.platform=="cygwin":
from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
from liblibra_core import *
import util.libutil as comn
def takeFirst(elem):
return elem[0]
def takeSecond(elem):
return elem[1]
[docs]def compute_oscillator_strengths(Hvib, Hprime_x, Hprime_y, Hprime_z, params):
"""Auxiliary function for the computation of the absorption/emission spectrum of the system
https://en.wikipedia.org/wiki/Oscillator_strength
Args:
Hvib ( list of CMATRIX ): "vibronic" Hamiltonians along the trajectory [units: a.u.]
Hprime_x ( list of CMATRIX ): i*hbar*<i|p_x|j> matrices along the trajectory [units: a.u.]
Hprime_y ( list of CMATRIX ): i*hbar*<i|p_y|j> matrices along the trajectory [units: a.u.]
Hprime_z ( list of CMATRIX ): i*hbar*<i|p_z|j> matrices along the trajectory [units: a.u.]
Returns:
[ [dE, f_ij, istep, i, j] ...], where:
* dE ( double ): energy gap [units: a.u.]
* f_ij ( double ): oscillator strength for a given transition [units: a.u.]
* istep ( int ): index of the time step along the given trajectory
* i ( int ): index of the initial (from where) state
* j ( int ): index of the final (to where) state
"""
homo = params["HOMO"]
nsteps = len(Hvib)
nmo = Hvib[0].num_of_cols
res = []
for istep in range(0,nsteps): # All data points
for i in range(0,homo+1): # All occupied orbitals
for j in range(homo+1, nmo): # All unoccupied orbitals
dE = Hvib[istep].get(j,j).real - Hvib[istep].get(i,i).real # >= 0.0
X = Hprime_x[istep].get(i,j)
X2 = (X.conjugate() * X ).real
Y = Hprime_y[istep].get(i,j)
Y2 = (Y.conjugate() * Y ).real
Z = Hprime_z[istep].get(i,j)
Z2 = (Z.conjugate() * Z ).real
f_ij = 0.0
if dE>0.0:
f_ij = (2.0/3.0)*( X2 + Y2 + Z2 )*(1.0/dE)
res.append([dE, f_ij, istep, i, j])
return res
[docs]def compute_spectrum(Hvib, Hprime_x, Hprime_y, Hprime_z, params):
"""Auxiliary function for the computation of the absorption/emission spectrum of the system
https://en.wikipedia.org/wiki/Oscillator_strength
Args:
Hvib ( list of CMATRIX ): "vibronic" Hamiltonians along the trajectory [units: a.u.]
Hprime_x ( list of CMATRIX ): i*hbar*<i|p_x|j> matrices along the trajectory [units: a.u.]
Hprime_y ( list of CMATRIX ): i*hbar*<i|p_y|j> matrices along the trajectory [units: a.u.]
Hprime_z ( list of CMATRIX ): i*hbar*<i|p_z|j> matrices along the trajectory [units: a.u.]
Returns:
tuple: ( E, osc_str, istep, init_st, fin_st ), where are the sorted data:
* E ( list of doubles ): energy gaps [units: a.u.]
* osc_str ( list of doubles ): oscillator strengths [units: a.u.]
* istep ( list of ints ): indices of the time steps
* init_st ( int ): indices of the initial states
* fin_st ( int ): indices of the final states
"""
opt = params["opt"]
res = compute_oscillator_strengths(Hvib, Hprime_x, Hprime_y, Hprime_z, params)
# sort list with key
sorted_res = None
if opt==0: # sort by energies
sorted_res = sorted(res, key=takeFirst)
elif opt==1: # sort by the magnitude of oscillator strength
sorted_res = sorted(res, key=takeSecond)
# Unpack the data into separate lists
E, osc_str, istep, init_st, fin_st = [], [], [], [], []
for it in sorted_res:
E.append(it[0])
osc_str.append(it[1])
istep.append(it[2])
init_st.append(it[3])
fin_st.append(it[4])
return E, osc_str, istep, init_st, fin_st
[docs]def get_step2_mb_sp_properties( params ):
"""
This function extracts information from es output files. It currently works for
cp2k, dftb+, gaussian
Args:
params ( dictionary ): parameters controlling the function execution
* **params["ks_orbital_indicies"] (list) - ks orbitals that were used in step2. Ex) [15,16,17, ... , 32,33,34]
* **params["logfile_directory"]** (string) - "../excitation_analysis/all_logfiles"
* **params["es_software"]** (string) - "cp2k" or "dftb+", etc
* **params["isUKS"]** (int) - 0 for spin-unpolarized, 1 for spin-polarized
* **params["number_of_states"]** (int) - number of ci states to extract
* **params["tolerance"]** (float) - cutoff for SD contribution
* **params["start_time"]** (int) - file index to start reading from
* **params["finish_time"]** (int) - file index to stop reading from
Returns:
Overlaps of SD at each step: S_sd_job
a time series of matrices
Time-overlaps of SD at each step: St_sd_job
a time series of Hvib matrices
One of each of the SP transitions (and its spin) that made up the considered CI states: sd_basis_states_unique
ex) [ ['1, 2'], 'alp' ] # 1 = homo, 2 = lumo
SP transitions for each excitation for each step: ci_basis_states_job
ex) [ ['1 ,2'] ]
CI coefficients for each excitation for each step: ci_coefficients_job
[ [1.0] ]
Spin components (alpha or beta excitaiton?) for each excitation for each step: spin_components_job
[ ['alp'] ]
"""
critical_params = ["logfile_directory", "es_software", "number_of_states", "start_time", "finish_time"]
default_params = {"isUKS":0, "tolerance":0.01 }
comn.check_input(params, default_params, critical_params)
start_time = params["start_time"]
finish_time = params["finish_time"]
curr_step = start_time
# Define ks_orbital_indicies based on the min_band and max_band
# Update params with ks_orbital_indicies
params["curr_step"] = curr_step
params["logfile_name"] = "step_"+str(params["curr_step"])+".out"
S_sd_job, St_sd_job = [], []
ci_energies_job, ci_basis_states_job, ci_coefficients_job = [], [], []
spin_components_job = []
sd_basis_states_unique = []
excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = step2_many_body.get_excitation_analysis_output( params )
ci_coefficients_raw_norm = step2_many_body.normalize_ci_coefficients(ci_coefficients_raw_unnorm)
print( excitation_energies )
print( ci_basis_raw )
print( ci_coefficients_raw_unnorm )
print( spin_components )
print( ci_coefficients_raw_norm )
# Now append the extracted excitation analysis output to the respective lists
ci_basis_states_job.append( ci_basis_raw )
ci_coefficients_job.append( ci_coefficients_raw_norm )
ci_energies_job.append( excitation_energies )
spin_components_job.append( spin_components )
# Extract the uniquie SD basis states from the ci basis states
for ci_basis_state_index in range( len( ci_basis_raw ) ):
for sd_basis_state_index in range( len( ci_basis_raw[ ci_basis_state_index ] ) ):
sd_basis_state_and_spin = [ ci_basis_raw[ ci_basis_state_index ][ sd_basis_state_index ] , spin_components[ ci_basis_state_index ][ sd_basis_state_index ] ]
if sd_basis_state_and_spin not in sd_basis_states_unique:
sd_basis_states_unique.append( sd_basis_state_and_spin )
#print( "Slater determinant basis states = ", sd_basis_states_unique )
curr_step += 1
# All other steps after initial step for this job
for step in range( finish_time-start_time-1 ):
params.update({ "curr_step":curr_step })
excitation_energies, ci_basis_raw, ci_coefficients_raw_unnorm, spin_components = step2_many_body.get_excitation_analysis_output( params )
ci_coefficients_raw_norm = step2_many_body.normalize_ci_coefficients(ci_coefficients_raw_unnorm)
# Extract the uniquie SD basis states from the ci basis states
for ci_basis_state_index in range( len( ci_basis_raw ) ):
for sd_basis_state_index in range( len( ci_basis_raw[ ci_basis_state_index ] ) ):
sd_basis_state_and_spin = [ ci_basis_raw[ ci_basis_state_index ][ sd_basis_state_index ] , spin_components[ ci_basis_state_index ][ sd_basis_state_index ] ]
if sd_basis_state_and_spin not in sd_basis_states_unique:
sd_basis_states_unique.append( sd_basis_state_and_spin )
#print( "Slater determinant basis states = ", sd_basis_states_unique )
# Now append the extracted excitation analysis output in _job variables
ci_basis_states_job.append( ci_basis_raw )
ci_coefficients_job.append( ci_coefficients_raw_norm )
ci_energies_job.append( excitation_energies )
spin_components_job.append( spin_components )
curr_step += 1
return S_sd_job, St_sd_job, sd_basis_states_unique, ci_basis_states_job, ci_coefficients_job, ci_energies_job, spin_components_job