# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

# # Rutile TiO$_2$ polaron hopping
# ## Overview
# In this notebook, we will:
# - Sample the potential energy surface of the polaron hop of TiO$_2$
# - Fit it to harmonic wells with both a hard cut off to exclude anharmonic points, or Boltzmann weighting of the fit
# - Calculate and compare the polaron hopping rate from Marcus theory

using CarrierCapture
using Plots, LaTeXStrings
using DataFrames
@show VERSION

# ## Adiabatic PES
# Here, we read the energies of linearly interpolation geometries from a set of DFT calculations. We then normalise it, such that its minimum sits at 0.

# +
pot = pot_from_file("data/tio2_interpolation.dat")

pot.QE_data.E = pot.QE_data.E .- minimum(pot.QE_data.E)

pot.Q0 = 0
pot.E0 = 0
# -

# Now, we fit the surface to a spline, by default of order 2.

pot.func_type = "spline"
fit_pot!(pot)

# ## Diabatic PESs
# We now split the potential in half. In this case, the saddle point is included, therefore we discard the maximum, which does not belong to the initial or final sets of sample points.

# +
pot_i, pot_f = cleave_pot(pot,discard_max=true)

# option to filter out points above kT
#filter_sample_points!(pot_i)
#filter_sample_points!(pot_f)

pot_i.T_weight = true
pot_f.T_weight = true

fit_pot!(pot_i)
fit_pot!(pot_f)
# -

# ## Plot
# Let's put it all together in one plot

plt = plot(pot.Q,pot.E,lab = "Adiabatic PES",thickness_scaling = 1.5)
plt = plot!(pot_i.Q,pot_i.E,lab = "Diabatic PES i")
plt = plot!(pot_f.Q,pot_f.E,lab = "Diabatic PES f")
scatter!(pot.QE_data.Q,pot.QE_data.E, lab="DFT Data",
    xlabel = "Fractional reaction coordinate",
    ylabel = "Relative energy (eV)",
    legend = true,
    title = "TiO2 <NAME>")


# ## Get transport properties

# +
tc = TransferCoord(pot_i, pot_f, pot)
coupling = get_coupling(tc)
reorg_ener = get_reorg_energy(tc)
activation_ener = get_activation_energy(tc)
transfer_rate = get_transfer_rate(tc)
dist = 2.904e-8
temp = 293.15
mobility = einstein_mobility(transfer_rate, 1, dist, temp)

println("Coupling: ", coupling, " eV")
println("Reorganisation enrgy: ", reorg_ener, " eV")
println("Activation energy: ", activation_ener, " eV")
println("Transfer rate: ", transfer_rate, " 1/s")
println("Einstein mobility: ", mobility, " cm^2V")
# -

# ## Comparison with NEB
# We repeat the whole procedure in one go for NEB results, and compare.

# +
# read in data
pot_neb = pot_from_file("data/tio2_neb.dat")

# normalise
pot_neb.QE_data.E = pot_neb.QE_data.E .- minimum(pot_neb.QE_data.E)

pot_neb.Q0 = 0
pot_neb.E0 = 0

# fit adiabatic surface
pot_neb.func_type = "spline"
fit_pot!(pot_neb)

# get diabatic surfaces
pot_i_neb, pot_f_neb = cleave_pot(pot_neb,discard_max=true)


pot_i_neb.T_weight = true
pot_f_neb.T_weight = true

fit_pot!(pot_i_neb)
fit_pot!(pot_f_neb)

# extract properties
tc = TransferCoord(pot_i_neb, pot_f_neb, pot_neb)
coupling = get_coupling(tc)
reorg_ener = get_reorg_energy(tc)
activation_ener = get_activation_energy(tc)
transfer_rate = get_transfer_rate(tc)
dist = 2.904e-8
temp = 293.15
mobility = einstein_mobility(transfer_rate, 1, dist, temp)

println("Coupling: ", coupling, " eV")
println("Reorganisation enrgy: ", reorg_ener, " eV")
println("Activation energy: ", activation_ener, " eV")
println("Transfer rate: ", transfer_rate, " 1/s")
println("Einstein mobility: ", mobility, " cm^2V")
# -

plt = plot(pot.Q,pot.E,lab = "Adiabatic PES interpolation",linewidth = 3,linecolor=1)
plt = plot!(pot_i.Q,pot_i.E,lab = "Diabatic PES interpolation",linewidth = 2,linecolor=1)
plt = plot!(pot_f.Q,pot_f.E,lab="",linewidth = 2,linecolor=1)
scatter!(pot.QE_data.Q,pot.QE_data.E, lab="DFT Data interpolation",
    markercolor=1,
    xlabel = "Fractional reaction coordinate",
    ylabel = "Relative energy (eV)",
    legend = true,
    title = "<NAME>")
plt = plot!(pot_neb.Q,pot_neb.E,lab = "Adiabatic PES NEB",linewidth = 3,linecolor=2,thickness_scaling = 1)
plt = plot!(pot_i_neb.Q,pot_i_neb.E,lab = "Diabatic PES NEB",linewidth = 2,linecolor = 2)
plt = plot!(pot_f_neb.Q,pot_f_neb.E,lab="",linewidth = 2,linecolor = 2)
scatter!(pot_neb.QE_data.Q,pot_neb.QE_data.E, lab="DFT Data NEB",
    markercolor=2,
    xlabel = "Fractional reaction coordinate",
    ylabel = "Relative energy (eV)",
    legend = true,
    title = "<NAME>")

StarcoderdataJupyterScriptsDedupFiltered