Tutorial on the Reproduction of selected Articles

%load_ext autoreload
%autoreload 2

# Save flag: Set to True to enable saving results (currently unused in this script)
save = False

# Verbose flag: Set to True to enable detailed logging
verbose = False

Setup

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Installation of QuantumDNA
# --------------------------

from importlib.util import find_spec

qDNA_installed = find_spec('qDNA') is not None

if not qDNA_installed:
    %pip install qDNA
    print("Successfully installed the 'qDNA' package.")
else:
    print("Package 'qDNA' is already installed.")

if verbose:
    %pip show qDNA

# ------------------------
# Directory Setup
# ------------------------

import os

# Use the current working directory as the root
ROOT_DIR = os.getcwd()

# Define directories to load data
DATA_DIR = os.path.join(ROOT_DIR, "data", "my_geometries")
os.makedirs(DATA_DIR, exist_ok=True)

# Define directory to save figures
SAVE_DIR = os.path.join(ROOT_DIR, "my_figures")
os.makedirs(SAVE_DIR, exist_ok=True)

if verbose:
    print(f"Data directory: '{DATA_DIR}' is ready.")
    print(f"Save directory: '{SAVE_DIR}' is ready.")

Package 'qDNA' is already installed.

Giese 1999 and Simserides 2014

[1] B. Giese et al., Angew. Chem. Int. Ed. 38, 996 (1999)

[2] C. Simserides, Chemical Physics 440, 31 (2014)

# Giese experiment: eta = 1.7, beta = 0.07

from qDNA import TB_Ham, DNA_Seq

tb_model_name = 'WM'
kwargs = dict(description='1P', particles=['hole'], unit='rad/fs', relaxation=False, source='Simserides2014')

GGG_per_G = []
for num_steps in np.arange(1, 9):
    bridge = 'TT' + 'GTT' * (num_steps-1)
    upper_strand = 'ACGCACGTCGCATAATATTAC' + 'G' + bridge + 'GGG' + 'TATTATATTACGC'
    tb_ham = TB_Ham(DNA_Seq(upper_strand, tb_model_name), **kwargs)

    donor_site = '(0, 21)' # guanine left from the bridge
    acceptor_sites = [f'(0, {21+len(bridge)+1})', f'(0, {21+len(bridge)+2})',  f'(0, {21+len(bridge)+3})']

    donor_avg_pop = tb_ham.get_average_pop(donor_site, donor_site)['hole']
    acceptor_avg_pop = 0
    for acceptor_site in acceptor_sites:
        acceptor_avg_pop += tb_ham.get_average_pop(donor_site, acceptor_site)['hole']

    GGG_per_G_ratio = acceptor_avg_pop / donor_avg_pop
    GGG_per_G.append(GGG_per_G_ratio)

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6.4, 4.8))

colors = [(0.2980392156862745, 0.4470588235294118, 0.6901960784313725),
 (0.3333333333333333, 0.6588235294117647, 0.40784313725490196)]

ax.plot(np.arange(1, 9) * 3 * 3.4, np.array(GGG_per_G)/(2*np.pi), 'o-', color=colors[0])

top_ticks = np.arange(0,100,10)
top_labels = np.arange(0,100,10)
ax.set_xticks(top_ticks)
ax.set_xticklabels(top_labels)

ax.set_yscale('log')
ax.set_xlim(0, 9*3*3.4)
ax.set_xlabel(r'distance [$\AA$]')
ax.set_ylabel('GGGperG23')

ax2 = ax.twiny()
# ax2.plot(np.arange(1, 9), np.array(ratio_list)/(2*np.pi) * np.e**2, 'o-', color=colors[1] )

ax2.set_xlabel(r'TT steps')
top_ticks = np.arange(0, 10)
top_labels = np.arange(0, 10)
ax2.set_xticks(top_ticks)
ax2.set_xticklabels(top_labels)

ax.set_yticks([])
ax2.set_yticks([])
ax.set_yticklabels([])
ax2.set_yticklabels([])

if save:
    fig_filename = input("Filename for Saving: ")
    plt.savefig(os.path.join(SAVE_DIR, fig_filename + '.pdf'))

Giese 2001 and Simserides 2014

[1] B. Giese, Nature 412, 318 (2001)

[2] C. Simserides, Chemical Physics 440, 31 (2014)

# experiment: beta = 0.6
# theory: beta = 0.8, beta = 0.07

from qDNA import TB_Ham, DNA_Seq, get_pop_fourier

def calc_difference(upper_strand, t):
    """ t in fs. """
    tb_ham = TB_Ham(DNA_Seq(upper_strand, 'WM'), description='1P', particles=['hole'], unit='rad/fs', relaxation=False, source='Simserides2014')
    difference = 0
    acceptor_sites = [f'(0, {len(bridge)+1})', f'(0, {len(bridge)+2})', f'(0, {len(bridge)+3})']
    for acceptor_site in acceptor_sites:
        amplitudes_dict, frequencies_dict, average_pop_dict = tb_ham.get_fourier('(0, 0)', acceptor_site, ["amplitude", "frequency", "average_pop"])
        amplitudes, frequencies, average_pop = amplitudes_dict['hole'], frequencies_dict['hole'], average_pop_dict['hole']
        difference += average_pop - get_pop_fourier(t, average_pop, amplitudes, frequencies)
    return average_pop, difference


def calc_t_mean(upper_strand, t_max, t_min=0, t_step=1):
    for t in range(t_min, t_max, t_step):
        average_pop, difference = calc_difference(upper_strand, t)
        if difference <= 0:
            return average_pop, t
    return average_pop, f"Mean population not reached in {t_max} fs"

t_mean_list, average_pop_list, transfer_rate_list = [], [], []
t_min_list = [2000, 32000, 225000, 250000, 140000, 3000, 0, 1000, 4000, 1000, 0, 1000, 4000, 3000, 3000, 1000]
t_steps_list = [1, 10, 10, 10, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

for i, num_T in enumerate(np.arange(1,17)):
    bridge = 'T' * num_T
    upper_strand = 'G' + bridge + 'GGG' + 'TATTATATTACGC'
    average_pop, t_mean = calc_t_mean(upper_strand, 10_000_000, t_min=t_min_list[i], t_step=t_steps_list[i])
    average_pop_list.append(average_pop)
    t_mean_list.append(t_mean)
    transfer_rate_list.append(average_pop/t_mean)

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

ax.plot(np.arange(3, 19)*3.4, np.array(transfer_rate_list)/(2*np.pi), 'o-')
ax.set_yscale('log')
ax.set_xlim(0, 70)
ax.set_xlabel(r'd [$\AA$]')
ax.set_ylabel('k [PHz]')

if save:
    fig_filename = input("Filename for Saving: ")
    plt.savefig(os.path.join(SAVE_DIR, fig_filename + '.pdf'))

Bittner 2006 & 2007

[1] E. Bittner, The Journal of Chemical Physics 125, 094909 (2006)

[2] E. Bittner, Journal of Photochemistry and Photobiology A: Chemistry 190, 328 (2007)

from qDNA import get_me_solver, plot_pops_heatmap

upper_strand = 'AAAAAAAAAA'
tb_model_name = 'ELM'
kwargs = dict(source='Bittner2007', unit='eV', relaxation=False, coulomb_interaction=2.5,
              exchange_interaction=1, init_e_state='(0, 4)', init_h_state='(0, 4)', t_end=300, t_unit='fs', t_steps=300)

me_solver = get_me_solver(upper_strand, tb_model_name, **kwargs)

fig, ax = plot_pops_heatmap(me_solver, heatmap_type='seaborn', vmax_list=[0.25, 0.25, 0.1])

if save:
    fig_filename = input("Filename for Saving: ")
    plt.savefig(os.path.join(SAVE_DIR, fig_filename + '.pdf'))

fig, ax = plot_pops_heatmap(me_solver, heatmap_type='matplotlib', vmax_list=[0.25, 0.25, 0.1])

if save:
    fig_filename = input("Filename for Saving: ")
    plt.savefig(os.path.join(SAVE_DIR, fig_filename + '.pdf'))

Mantela 2023

[1] M. Mantela, Phys. Chem. Chem. Phys. 25, 7750 (2023)

from qDNA import load_xyz, Base, BasePair, Dimer

xyz_identifier, xyz_data = load_xyz("A1", os.path.join(DATA_DIR, "Mantela"))
base1 = Base(xyz_identifier, xyz_data)

xyz_identifier, xyz_data = load_xyz("c1'", os.path.join(DATA_DIR, "Mantela"))
base2 = Base(xyz_identifier, xyz_data)

xyz_identifier, xyz_data = load_xyz("A2", os.path.join(DATA_DIR, "Mantela"))
base3 = Base(xyz_identifier, xyz_data)

xyz_identifier, xyz_data = load_xyz("c2'", os.path.join(DATA_DIR, "Mantela"))
base4 = Base(xyz_identifier, xyz_data)

basepair1 = BasePair(base1, base2)
basepair2 = BasePair(base3, base4)

dimer = Dimer(basepair1, basepair2)

print("TB parameters")
print("-------------------------------")
print(f"E_HOMO1: {dimer.molecule1.E_HOMO}")
print(f"E_LUMO1: {dimer.molecule1.E_LUMO}")
print("-------------------------------")
print(f"E_HOMO2: {dimer.molecule2.E_HOMO}")
print(f"E_LUMO2: {dimer.molecule2.E_LUMO}")
print("-------------------------------")
print(f"t_HOMO: {dimer.t_HOMO}")
print(f"t_LUMO: {dimer.t_LUMO}")

TB parameters
-------------------------------
E_HOMO1: -8.429037928082415
E_LUMO1: -4.432891842321037
-------------------------------
E_HOMO2: -8.430148205528646
E_LUMO2: -4.43267272482506
-------------------------------
t_HOMO: -0.035940257425427
t_LUMO: -0.09015171853086595