5
$\begingroup$

I am thinking about calculating the HOMO-LUMO energy gap for $\ce{N2}$ molecule.

I thought I could model this as a particle in a box problem to find $E_\text{HOMO}$ and $E_\text{LUMO}$.

The formula I used was $E_n = \frac{h^2 n^2}{8 m_e L^2}$

I put in $L$ as the bond length of nitrogen. This gave me the wrong results. I hear this kind of thing works for conjugated molecules, and I thought it would work for $\ce{N2}$.

Is there a simple hack of this variety that can be used to calculate this or generate a short table of such transitions?

I think DFT is good for ground state-related affairs and besides needs a bit more setting up.

Is there also maybe a reference that does this kind of thing for $\ce{N2}$ molecule nicely with full detail? I want to use the differences to predict the wavelengths to be observed by UV spectroscopy by $\Delta E = h \frac{c}{\lambda}$.

I have to put my pride to the side and also trust that this is a nice good safe space for my code.

Here is what I have tried. I wanted to calculate the ground states, but it is not quite working. Also, I think the simple 1-d DFT only gets the ground state, so my scheme of doing E[1] - E[0] is wrong.

I probably need to upgrade this to TDDFT.

I want to plot transition probabilities vs wavelength at the end.

This code is written in Python. Can someone help with the missing pieces?

I have now attempted to generate excited states using ordinary perturbation theory from the ground states.

It generates a plot but it is not showing the trend I would expect. I expect the HOMO-LUMO Gap to fall, but it is not falling.

My code is wrong? Or my chemistry is horrible.

About the chemistry, I hacked my way to try to create an approximation of the TD DFT without time.

# Kevin Njokom
# Written completely from scratch


import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg  as linear_algebra
import math as math
import periodictable as elements




# Units
# hbar = e = m = 1


# constants
# number_of_electrons = 14  # Number of electrons
size = 50 #150  # Number of grid points (increased for accuracy) # at least 4
tolerance = 1e-40  # A small epsilon
max_iterations = 100



# Constants
h = 6.62607015e-34  # Planck's constant in J s
c = 299792458       # Speed of light in m/s
au_to_joules = 4.35974417e-18  # Conversion factor from atomic units to joules
angstrom_to_meter = 1e-10  # Conversion factor from angstroms to meters


# Global lists

delta_Energy_keeper =[]
psi_ground_keeper =[]
psi_excited_keeper = []

Vectoria_2 = []
Excited_energy_2 = []


def give_me_angstroms(energy_au):
    # Convert energy from atomic units to joules
    energy_joules = energy_au * au_to_joules 
    # Calculate wavelength in meters
    wavelength_meter = h * c / energy_joules
    # Convert wavelength to angstroms
    wavelength_angstrom = wavelength_meter / angstrom_to_meter
    # trying something
    return wavelength_angstrom

                   
# Finite difference approximation for the second derivative implemented in code

def create_second_derivative_operator(size):

    # Create the main diagonal
    main_diag = np.ones(size) * -2
    
    # Create the sub-diagonals
    sub_diag = np.ones(size-1)
    
    # Construct the tridiagonal matrix
    second_derivative_operator = np.diag(main_diag) + np.diag(sub_diag, k=1) + np.diag(sub_diag, k=-1)
    
    return second_derivative_operator
    
# Definition of the Reimann Integral implemented in code

def reimann_integration(x_co_ordinate, function_to_integrate):
    finite_difference =x_co_ordinate[1] - x_co_ordinate[0]
    Integrate = np.sum(function_to_integrate * finite_difference)
    Integrated = Integrate
    return Integrated

# This normalizes a vector

def vector_normalization(vector):
    magnitude = np.linalg.norm(vector)
    normalized_vector = vector / magnitude
    return normalized_vector

# poisson equation solver
# second_derivative_operator * potential  = electron density

def poison_equation_solver( second_derivative_operator, electron_density):
    potential = np.linalg.tensorsolve(second_derivative_operator,electron_density)
    return potential


# create wave function

def create_initial_wave_function(size,number_of_electrons):
    dim = size
    initial_wave_function = np.ones(dim)
    return initial_wave_function


def single_electron_wave_function(size):
    Hamiltonian = -0.5 *create_second_derivative_operator(size)
    eigen_values, eigen_vectors = linear_algebra.eigh(Hamiltonian)
    return eigen_vectors

def single_electron_energies(size):
    Hamiltonian = -0.5 *create_second_derivative_operator(size)
    eigen_values, eigen_vectors = linear_algebra.eigh(Hamiltonian)
    return eigen_values

# each wavefunction is in an orbit
# each contributes a density
# all the densities are added and multiplied by 2
# This is the electron density

def electron_density(vector_of_vectors):
    # Square each element of the vectors and sum them up
    sum_of_element_in_vector_of_square_vectors = np.array(list(map(lambda vec: np.sum(np.square(vec)), vector_of_vectors)))
    return 2*(sum_of_element_in_vector_of_square_vectors)

# Initial electron density guess

electron_density_guess = electron_density(single_electron_wave_function(size))

# Achievement unlocked! I have independently computed the hatree potential alone from scratch! No peeking, no reference!
# Hartree potential here is computed by pseudo -inverse methods

def hartree_potential(size, number_of_electrons):
    eigen_vectors = single_electron_wave_function(size)
    vector_of_normalized_vectors = np.array(list(map(lambda vec: vector_normalization(vec), eigen_vectors)))
    new_electron_density = electron_density(vector_of_normalized_vectors)
    Apinv = np.linalg.pinv(create_second_derivative_operator(size))  # Finds pseudo-inverse
    V_hartree_potential = Apinv.dot(new_electron_density*(np.ones(size)))
    V_hartree_potential = V_hartree_potential*number_of_electrons
    return V_hartree_potential

# Exchange potential

def exchange_potential(size, number_of_electrons):
    eigen_vectors = single_electron_wave_function(size)
    vector_of_normalized_vectors = np.array(list(map(lambda vec: vector_normalization(vec), eigen_vectors)))
    new_electron_density = electron_density(vector_of_normalized_vectors)   
    V_exchange_potential  = -(3. / np.pi * new_electron_density)**(1. / 3.)
    return V_exchange_potential


def complete_hamiltonian_with_interaction(size, number_of_electrons):
    kinetic_contribution = -0.5 *create_second_derivative_operator(size)
    complete_hamiltonian_with_interaction = -0.5 *kinetic_contribution + np.diagflat(hartree_potential(size,number_of_electrons)) + np.diagflat(exchange_potential(size,number_of_electrons))
    eigen_values,eigen_vectors = linear_algebra.eigh(complete_hamiltonian_with_interaction)
    return eigen_vectors

def complete_hamiltonian_with_interaction_energies(size, number_of_electrons):
    kinetic_contribution = -0.5 *create_second_derivative_operator(size)
    complete_hamiltonian_with_interaction = kinetic_contribution + np.diagflat(hartree_potential(size,number_of_electrons)) + np.diagflat(exchange_potential(size,number_of_electrons))
    eigen_values,eigen_vectors = linear_algebra.eigh(complete_hamiltonian_with_interaction)
    return eigen_values

def print_separator(character='-', length=50):
    separator = character * length
    print(separator)


def make_matrix(omega, size):
    # Create a matrix of zeros of size size by size
    mat = np.zeros((size, size))
    # Set diagonal elements to -omega
    np.fill_diagonal(mat, -omega)
    return mat
    

def matrix_form_solve_for_x_Ax_omega_x(vector1, vector2, vector3):
    # Perform the matrix multiplication vector1 . matrix
    Omega_number = np.dot(vector1[0], vector3[0])
    # Perform the matrix multiplication result1 . vector2^T
    Omega_number = Omega_number * vector2[0]
    return Omega_number  #Actually, this is a number


def self_consistency_loop(size, number_of_electrons, max_iterations):
    electron_density_guess = electron_density(single_electron_wave_function(size))
    nuevo_eigen_vectors = complete_hamiltonian_with_interaction(size, number_of_electrons)
    normalized_nuevo_eigen_vectors = vector_normalization(nuevo_eigen_vectors)
    nuevo_electron_density = electron_density(normalized_nuevo_eigen_vectors)
    old_electron_density = electron_density_guess
    Energetica = []
    Vectoria = []
    leading_energy_correction = []
    leading_energy_correction_2 = []
    iteration = 0

    while iteration < max_iterations:
        if abs(nuevo_electron_density[0] - old_electron_density[0]) < tolerance:
            print("new electron density" + " " + str(nuevo_electron_density[0]) + "\n")
            print("convergence achieved" + "\n")
            print("Energies" + str(Energetica) + "\n")
            print("Energy Difference" + " " + str(Energetica[1] - Energetica[0]) + "\n")
            print("tolerance to within " + str(abs(nuevo_electron_density[0] - old_electron_density[0])) + "\n")
            print("Iteration number " + str(iteration) + "\n")
            print_separator()
            leading_energy_correction = matrix_form_solve_for_x_Ax_omega_x(Vectoria,V_exchange_potential,Vectoria)
            leading_energy_correction_2.append(matrix_form_solve_for_x_Ax_omega_x(Vectoria,V_exchange_potential,Vectoria))
            sorted_Energies = sorted(Energetica)
            Excited_energy_n, Excited_state = linear_algebra.eigh(make_matrix(leading_energy_correction + sorted_Energies[0],size))
            delta_Energy = leading_energy_correction
            delta_Energy_keeper.append(delta_Energy)
            diff_now = delta_Energy_keeper[-1]
            wavelength_lambda = give_me_angstroms(diff_now)
            delta_Energy_keeper.append(delta_Energy)
            psi_ground_keeper.append(Vectoria[0])
            psi_excited_keeper.append(Excited_state[0])
            prob_trans = (np.dot(psi_ground_keeper[-1], psi_excited_keeper[-1]))**2
            wavelength_lambda = give_me_angstroms(abs(diff_now))
            return [wavelength_lambda, prob_trans, diff_now]

        # Update the electron density and energies if the convergence criterion is not met
        old_electron_density = nuevo_electron_density
        V_exchange_potential = -(3. / np.pi * old_electron_density) ** (1. / 3.)
        Apinv = np.linalg.pinv(create_second_derivative_operator(size))
        V_hartree_potential = Apinv.dot(old_electron_density * (np.ones(size)))
        V_hartree_potential = V_hartree_potential * number_of_electrons
        kinetic_contribution = -0.5 *create_second_derivative_operator(size)
        newestHamiltonian = kinetic_contribution + np.diagflat(V_hartree_potential) + np.diagflat(V_exchange_potential)
        new_Energies, nuevo_eigen_vectors = linear_algebra.eigh(newestHamiltonian)
        Energetica = new_Energies
        Vectoria = nuevo_eigen_vectors
        normalized_nuevo_eigen_vectors = vector_normalization(nuevo_eigen_vectors)
        nuevo_electron_density = electron_density(normalized_nuevo_eigen_vectors)

        print("new electron density" + " " + str(nuevo_electron_density[0]) + "\n")
        print("convergence not achieved" + "\n")
        print("Energies" + str(Energetica) + "\n")
        print("Iteration number " + str(iteration) + "\n")
        print("tolerance to within " + str(abs(nuevo_electron_density[0] - old_electron_density[0])) + "\n")
        print_separator()
        print(Vectoria)
        
        iteration += 1  # Increment iteration after convergence check and updates
        print("Energy  = "+ str(np.min(Energetica)) +" " + "for iteration number" + " " + str(iteration))
    
# Define dictionary mapping atomic numbers to element symbols
atomic_number_to_symbol = {1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N', 8: 'O', 9: 'F', 10: 'Ne',
                           11: 'Na', 12: 'Mg', 13: 'Al', 14: 'Si', 15: 'P', 16: 'S', 17: 'Cl', 18: 'Ar',
                           19: 'K', 20: 'Ca', 21: 'Sc', 22: 'Ti', 23: 'V', 24: 'Cr', 25: 'Mn', 26: 'Fe',
                           27: 'Co', 28: 'Ni', 29: 'Cu', 30: 'Zn', 31: 'Ga', 32: 'Ge', 33: 'As', 34: 'Se',
                           35: 'Br', 36: 'Kr', 37: 'Rb', 38: 'Sr', 39: 'Y', 40: 'Zr', 41: 'Nb', 42: 'Mo',
                           43: 'Tc', 44: 'Ru', 45: 'Rh', 46: 'Pd', 47: 'Ag', 48: 'Cd', 49: 'In', 50: 'Sn',
                           51: 'Sb', 52: 'Te', 53: 'I', 54: 'Xe', 55: 'Cs', 56: 'Ba', 57: 'La', 58: 'Ce',
                           59: 'Pr', 60: 'Nd', 61: 'Pm', 62: 'Sm', 63: 'Eu', 64: 'Gd', 65: 'Tb', 66: 'Dy',
                           67: 'Ho', 68: 'Er', 69: 'Tm', 70: 'Yb', 71: 'Lu', 72: 'Hf', 73: 'Ta', 74: 'W',
                           75: 'Re', 76: 'Os', 77: 'Ir', 78: 'Pt', 79: 'Au', 80: 'Hg', 81: 'Tl', 82: 'Pb',
                           83: 'Bi', 84: 'Po', 85: 'At', 86: 'Rn', 87: 'Fr', 88: 'Ra', 89: 'Ac', 90: 'Th',
                           91: 'Pa', 92: 'U', 93: 'Np', 94: 'Pu', 95: 'Am', 96: 'Cm', 97: 'Bk', 98: 'Cf',
                           99: 'Es', 100: 'Fm', 101: 'Md', 102: 'No', 103: 'Lr', 104: 'Rf', 105: 'Db',
                           106: 'Sg', 107: 'Bh', 108: 'Hs', 109: 'Mt', 110: 'Ds', 111: 'Rg', 112: 'Cn',
                           113: 'Nh', 114: 'Fl', 115: 'Mc', 116: 'Lv', 117: 'Ts', 118: 'Og'}

# Define lists to store wavelengths and probabilities
wavelengths = [[] for _ in range(20)]
probabilities = [[] for _ in range(20)]
Energy_differences = [[] for _ in range(20)]

# Loop through different numbers of electrons
for num_electrons in range(1, 21):
    # Call self_consistency_loop function
    result = self_consistency_loop(size, num_electrons, max_iterations)
    # Store the result
    wavelengths[num_electrons - 1].append(result[0])
    probabilities[num_electrons - 1].append(result[1])
    Energy_differences[num_electrons - 1].append(result[2])
    

# Plotting
plt.figure(figsize=(10, 6))
for num_electrons in range(1, 21):
    if num_electrons == 14:
        label = 'Si or N2'  # Label as carbon or nitrogen molecule for 14 electrons
    else:
        label = atomic_number_to_symbol.get(num_electrons, f"Element {num_electrons}")  # Get symbol from dictionary
    # Plot data points
    plt.plot(wavelengths[num_electrons-1], probabilities[num_electrons-1], marker='o', linestyle='', label=label)
    # Add labels for each point
    for x, y in zip(wavelengths[num_electrons-1], probabilities[num_electrons-1]):
        if num_electrons == 14:
            plt.text(x, y, label, fontsize=8, ha='right', va='bottom')
        else:
            plt.text(x, y, atomic_number_to_symbol[num_electrons], fontsize=8, ha='right', va='bottom')

plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Probability of Transition')
plt.title('Probability of Transition vs Wavelength for Different Numbers of Electrons (atoms or molecules)')
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=5)
plt.grid(True)
plt.show()

print("---------------------------------")
print(delta_Energy_keeper)
print("---------------------------------")
print(psi_ground_keeper[0])
print("---------------------------------")
print(psi_excited_keeper[0])
print("---------------------------------")
print(Excited_energy_2)
print("---------------------------------")
print(Energy_differences)
print("---------------------------------")
$\endgroup$
3
  • 1
    $\begingroup$ While it is an interesting idea I suspect a 1 dimensional box is too far from the reality to get anything better than an order of magnitude estimate. But let's soldier on. Could you edit the question to show exactly the maths you did and the result you obtained; the result you expect with a reference would be useful as well. "The wrong results" really is nowhere near enough, we need to see what you did. $\endgroup$
    – Ian Bush
    Feb 16 at 10:18
  • $\begingroup$ For future reference: for the body of questions, answers, and comments, chemistry.se offers to use mhchem as a comfortable method to add chemical equations and report numerical values (\pu{}) including a non-breakable space. It should not to be used in the title of a post. $\endgroup$
    – Buttonwood
    Feb 16 at 15:14
  • $\begingroup$ Hello @IanBush , I have now shown progress thus far. I really hope this could maybe happen $\endgroup$ Feb 26 at 7:57

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.