Skip to content
Snippets Groups Projects
Commit a3340c0b authored by moritz.buchhorn's avatar moritz.buchhorn
Browse files

Wrap up week 7

parent 5b8c3a98
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
07-tddft/07-tddft.ipynb
+ 99
− 218
View file @ a3340c0b
%% Cell type:markdown id: tags:
# B.CTCT: Time-Dependent Density Functional Theory
Angeregte Zustände, Näherung, Intensität und Shift
Die überwiegende Mehrheit quantenchemischer Methoden beschäftigt sich mit dem elektronischen Grundzustand. Angeregte Zustände sind aber für die Berechnung von elektronischen Spektren und die Interpretation photochemischer Prozesse außerordentlich wichtig. Eine verbreitete und relativ verlässliche Methode für die Vorhersage von elektronischen Übergängen ist TDDFT. Mittels TDDFT können die Energien und Intensitäten elektronische Übergänge berechnet werden, die Rechnungen basieren dabei auf der Elektronendichte des Grundzustands. Es werden also nicht wirklich angeregten Zustände berechnet.
%% Cell type:markdown id: tags:
Wie häufig bei der Verwendung von Dichtefunktionaltheorie ist die quantitative Vorhersage elektronischer Übergänge schwierig. Es ist üblich, dass die Übergangsenergien verschoben sind; diese Verschiebung ist aber innerhalb einer Rechnung in der Regel systematisch. Ebenfalls können wie üblich Trends untersucht werden, d.h. ähnliche Systeme liefern ähnlich gute Ergebnisse.
Wir werden uns mit einer einfachen Reihe beschäftigen: lineare Polyene mit n Doppelbindungen, mit n von 2 bis 5. In der Übung beigefügt sind optimierte Strukturen der Verbindungen, sodass Sie direkt die TDDFT-Rechnungen ausführen können.
%% Cell type:code id: tags:
``` python
import plotly.graph_objects as go
import subprocess as sub
import support
```
%% Cell type:code id: tags:
``` python
vol = support.Volumetric.from_gbw_file("benzene-tddft.gbw", 3, "mo")
xyz = support.XYZ.from_xyz_file("polyene-3-opt.xyz")
xyz.get_molecular_viewer().show()
```
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` python
threshold = 6*10**-3
xyz = support.XYZ.from_xyz_file("benzene-opt.xyz")
xyz.get_molecular_viewer()
xyz.get_mesh_from_gbw("benzene-tddft.gbw", 0, mode="mo", threshold=threshold)
xyz.get_mesh_from_gbw("benzene-tddft.gbw", 0, mode="mo", threshold=-threshold)
# xyz.get_mesh_from_cube("benzene-tddft.mo5a.cube", threshold)
# xyz.get_mesh_from_cube("benzene-tddft.mo5a.cube", -threshold)
xyz.viewer.show()
```
Die einzelnen Rechnungen per Hand vorzubereiten und auszuführen wäre eine unnötig aufwändige und fehleranfällige Aufgabe, verwenden Sie deshalb die Schleife in der nächsten Zelle, um die _input files_ zu schreiben.
📝 Führen Sie die nächste Zelle aus und prüfen Sie in ihrem Ordner, ob die _input files_ korrekt erstellt wurden. Führen Sie dann die übernächste Zelle aus, um die Rechnungen zu starten.
%% Cell type:code id: tags:
``` python
def get_tddft_spectrum(output_file: str) -> tuple[list[int], list[float], list[float]]:
indices = []
wavenumbers = []
wavelengths = []
foscs = []
with open(output_file) as f:
for line in f:
if "ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS" in line:
next(f)
next(f)
next(f)
next(f)
while True:
line = next(f)
if len(line.split()) == 0: break
line = line.split()
indices.append(int(line[0]))
wavenumbers.append(float(line[1]))
wavelengths.append(float(line[2]))
foscs.append(float(line[3]))
return indices, wavenumbers, wavelengths, foscs
ns = [2,3,4,5]
for n in ns:
with open(f"polyene-{n}-tddft.inp", "w") as f:
f.write(f"""
! UKS BP def2-SVP
indices, wavenumbers, wavelengths, foscs = get_tddft_spectrum("benzene-tddft.out")
fig = go.Figure(layout=go.Layout(
title="Benzene TDDFT Spectrum",
# xaxis_title="Wavenumber / cm<sup>-1</sup>",
xaxis_title="Wavelength / nm",
yaxis_title="Oscillator Strength"))
# fig.add_trace(go.Bar(x=wavelengths, y=foscs, text=indices))
fig.add_trace(go.Scatter(x=wavelengths, y=foscs, text=indices, mode='text'))
fig.show()
%tddft nroots 25 end
* xyzfile 0 1 polyene-{n}-opt.xyz
""")
```
%% Cell type:code id: tags:
``` python
for n in ns:
print(f"Running polyene-{n}-tddft.inp")
sub.run(f"orca polyene-{n}-tddft.inp > polyene-{n}-tddft.out", shell=True)
```
from typing import Iterable
import numpy as np
%% Cell type:markdown id: tags:
class Volumetric:
def __init__(self,
origin: Iterable[float],
res_vec1: int, res_vec2: int, res_vec3: int,
vec1: Iterable[float],
vec2: Iterable[float],
vec3: Iterable[float],
volumetrics: Iterable[float]) -> None:
self.origin = np.array(origin)
self.res_vec1 = res_vec1
self.res_vec2 = res_vec2
self.res_vec3 = res_vec3
self.vec1 = np.array(vec1)
self.vec2 = np.array(vec2)
self.vec3 = np.array(vec3)
self.volumetrics = volumetrics
return
@classmethod
def read_cube_file(cls, file: str):
with open(file, "r") as f:
next(f)
next(f)
natoms, x0, y0, z0 = f.readline().split()
res_vec1, x1, y1, z1 = f.readline().split()
res_vec2, x2, y2, z2 = f.readline().split()
res_vec3, x3, y3, z3 = f.readline().split()
natoms, x0, y0, z0 = (abs(int(natoms)), float(x0), float(y0), float(z0))
res_vec1, x1, y1, z1 = (int(res_vec1), float(x1), float(y1), float(z1))
res_vec2, x2, y2, z2 = (int(res_vec2), float(x2), float(y2), float(z2))
res_vec3, x3, y3, z3 = (int(res_vec3), float(x3), float(y3), float(z3))
# skip all the atom coordinates
for i in range(natoms):
next(f)
if ".mo" in file:
next(f)
volumetrics = f.readlines()
volumetrics = [float(i) for line in volumetrics for i in line.split()]
return cls(
(x0, y0, z0),
res_vec1, res_vec2, res_vec3,
(x1, y1, z1),
(x2, y2, z2),
(x3, y3, z3),
volumetrics
)
def voxel_generator(self, lower: float, upper: float):
for x in range(self.res_vec1):
for y in range(self.res_vec2):
for z in range(self.res_vec3):
voxel = self.volumetrics[z + y * self.res_vec3 + x * (self.res_vec3 * self.res_vec2)]
if voxel > lower and voxel < upper:
coord = x * self.vec1 + y * self.vec2 + z * self.vec3
yield voxel, *coord
lower = 6*10**-5
vol = Volumetric.read_cube_file("benzene-tddft.mo5a.cube")
voxels = [voxel for voxel in vol.voxel_generator(lower, 1)]
voxels = np.array(voxels)
print(len(voxels))
print(voxels[:10])
```
&#x1f4dd; Öffnen Sie eines der _output files_. Gegen Ende der Datei finden Sie die berechneten elektronischen Übergänge mit ihrer Intensität (_fosc_, Oszillatorstärke):
<pre style="color:#0F0F0F; background-color:#F0F0F0; padding: 1em;">
-----------------------------------------------------------------------------
ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS
-----------------------------------------------------------------------------
State Energy Wavelength fosc T2 TX TY TZ
(cm-1) (nm) (au**2) (au) (au) (au)
-----------------------------------------------------------------------------
1 32438.4 308.3 0.000000000 0.00000 0.00000 0.00000 0.00000
2 36047.7 277.4 0.000000000 0.00000 0.00000 0.00000 0.00000
3 36550.0 273.6 0.000000000 0.00000 0.00000 -0.00000 0.00000
</pre>
Ein Stück darüber sind die Orbitale gezeigt, zwischen denen die jeweiligen Übergänge stattfinden:
<pre style="color:#0F0F0F; background-color:#F0F0F0; padding: 1em;">
STATE 1: E= 0.147800 au 4.022 eV 32438.4 cm**-1 <S**2> = 2.000000
23a -> 25a : 0.176683 (c= 0.42033655)
24a -> 26a : 0.316259 (c= 0.56236953)
23b -> 25b : 0.176683 (c= -0.42033655)
24b -> 26b : 0.316259 (c= -0.56236953)
</pre>
%% Cell type:code id: tags:
Welche Übergänge sind die intensivsten? Welche Orbitale sind beteiligt?
``` python
import sklearn.cluster as skc
%% Cell type:markdown id: tags:
cluster = skc.DBSCAN(eps=0.5).fit(voxels[:,1:4])
print(cluster.labels_)
print(len(cluster.labels_))
_ = [print(l, end=" ") for l in cluster.labels_]
print()
print(cluster.labels_.max())
```
&#x1f4dd; Plotten Sie das Linienspektrum mithilfe der nächsten Zelle. Sie müssen nichts ändern, die Daten werden direkt aus den _output files_ entnommen. Entspricht der gefundene Trend ihren Erwartungen?
%% Cell type:code id: tags:
``` python
i = 0
fig = go.Figure()
fig.add_trace(go.Mesh3d(x=voxels[cluster.labels_ == i][:,1],
y=voxels[cluster.labels_ == i][:,2],
z=voxels[cluster.labels_ == i][:,3],
opacity=0.5,
alphahull=1
))
fig.show()
```
fig = go.Figure(layout=go.Layout(
title="TDDFT Spectra",
xaxis_title="Wavelength / nm",
yaxis_title="Oscillator Strength"))
%% Cell type:code id: tags:
for n in ns:
spectrum = support.Spectrum.from_tddft_output(f"polyene-{n}-tddft.out")
wavelengths, intensities = spectrum.get_line_spectrum("nm")
# Fügt die Nummern hinzu
fig.add_trace(go.Scatter(
x=spectrum.get_energies_as("nm"),
y=spectrum.intensities,
text=spectrum.indices,
mode='text',
name=f"{n}-Polyen"
))
# Fügt die Linien hinzu
fig.add_trace(go.Scatter(
x=wavelengths,
y=intensities,
mode='lines',
name=f"{n}-Polyen"
))
``` python
fig = go.Figure()
for i in range(9):
cluster_i = voxels[cluster.labels_ == i]
fig.add_trace(go.Scatter3d(x=cluster_i[:,1],
y=cluster_i[:,2],
z=cluster_i[:,3],
mode='markers',
opacity=0.5
))
fig.show()
```
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` python
import pyvista as pv
Für die Deutung photochemischer Prozesse ist es sehr wichtig, herauszufinden, zwischen welchen Orbitalen elektronische Übergänge stattfinden. Da von niemandem erwartet werden kann, Orbitalkoeffizienten mental in eine elektronendichte zu übersetzen, werden Orbitale dreidimensional dargestellt.
xyz = np.ndarray((0, 3))
ijk = np.ndarray((0, 3))
fig = go.Figure()
for i in range(9):
cloud = pv.PolyData(voxels[cluster.labels_ == i][:,1:4])
surf = cloud.delaunay_3d(alpha=0.4, tol=0.001, offset=2.5)
surf = surf.extract_geometry()
print(surf.volume)
print(surf.volume)
surf = surf.triangulate()
print(len(surf.points))
print(surf.is_all_triangles)
surf = surf.decimate(0.9)
print(len(surf.points))
ijk = np.concatenate((ijk, surf.faces.reshape(-1, 4)[:, 1:] + len(xyz)))
xyz = np.concatenate((xyz, surf.points))
# xyz = surf.points
# ijk = surf.faces.reshape(-1, 4)[:, 1:]
fig.add_trace(
go.Mesh3d(x=xyz[:,0],
y=xyz[:,1],
z=xyz[:,2],
i=ijk[:,0],
j=ijk[:,1],
k=ijk[:,2],
opacity=0.5,
# alphahull=0
)
)
fig.show()
```
&#x1f4dd; Sehen Sie sich für ein beliebiges Molekül aus der Reihe den intensivsten Übergang im _output file_ an. Welche Orbitale sind beteiligt? Geben Sie die gefundenen Orbitalnummern in die Liste `mos` in der nächsten Zelle ein und starten Sie sie.
Hinweis: der Wert unter `threshold` gibt an, bei welchem Funktionswert die Orbitaloberfläche abgeschnitten wird. Je nach Orbital müssen Sie diesen Wert etwas anpassen.
%% Cell type:code id: tags:
``` python
import pyvista as pv
cloud = pv.PolyData(voxels[:,1:4])
surf = cloud.delaunay_3d(alpha=0.4, tol=0.001, offset=2.5)
# surf = cloud.reconstruct_surface(nbr_sz=40)
# surf = surf.extract_surface()
surf = surf.extract_geometry()
print(surf.volume)
# surf = surf.smooth_taubin(n_iter=100, edge_angle=5, feature_angle=10)
print(surf.volume)
surf = surf.triangulate()
print(len(surf.points))
print(surf.is_all_triangles)
surf = surf.decimate(0.9)
print(len(surf.points))
xyz = surf.points
ijk = surf.faces.reshape(-1, 4)[:, 1:]
# _ = [print(x) for i,x in enumerate(surf.faces) if i % 4 == 0]
surf.plot()
mos = [20, 21]
threshold = 5*10**-2
for mo in mos:
xyz = support.XYZ.from_xyz_file("polyene-4-opt.xyz")
xyz.make_mesh_from_gbw("polyene-4-tddft.gbw", mo, mode="mo", threshold=threshold, verbose=False, dynamic_threshold=True)
xyz.make_mesh_from_gbw("polyene-4-tddft.gbw", mo, mode="mo", threshold=-threshold, verbose=False, dynamic_threshold=True)
xyz.viewer.update_layout(title_text=f"Rendered molecular orbital {mo}")
xyz.viewer.show()
```
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` python
fig = go.Figure()
fig.add_trace(go.Mesh3d(x=xyz[:,0],
y=xyz[:,1],
z=xyz[:,2],
i=ijk[:,0],
j=ijk[:,1],
k=ijk[:,2],
opacity=0.5,
# alphahull=0
))
fig.show()
```
Eine kompakte Art, die unterschiedliche Dichte in Grund- und angeregtem Zustand darzustellen ist über die Differenzdichte. Hier werden die Elektronendichten von Grnd- und angeregtem Zustand voneinander abgezogen, sodass in einer dreidimensionalen Darstellung sichtbar ist, "von wo nach wo" der elektronische Übergang erfolgt. Üblicherweise wird der negative Teil (also der Elektronendichteverlust) rot dargestellt, und der positive blau.
&#x1f4dd; Suchen Sie sich im Spektrum der Reihe einen Übergang aus und benutzen Sie die nächste Zelle, um die dazugehörige Differenzdichte als 3D-Darstellung zu erhalten. Geben Sie die Übergänge, die Sie visualisieren wollen, in der Variable `transitions` ein.
Das Erstellen einer Differenzdichte kann etwas Zeit in Anspruch nehmen.
%% Cell type:code id: tags:
``` python
transitions = [23]
threshold = 5*10**-4
for transition in transitions:
xyz = support.XYZ.from_xyz_file("polyene-3-opt.xyz")
xyz.make_mesh_from_gbw("polyene-3-tddft.gbw", transition, mode="diffdens", threshold=threshold, verbose=True, dynamic_threshold=True)
xyz.make_mesh_from_gbw("polyene-3-tddft.gbw", transition, mode="diffdens", threshold=-threshold, verbose=True, dynamic_threshold=True)
xyz.viewer.show()
```
......
07-tddft/make_exercise.bash 0 → 100644
+ 6
− 0
View file @ a3340c0b
mkdir 07-tddft
cp 07-tddft.ipynb 07-tddft/
cp *.py 07-tddft/
cp polyene*opt.xyz 07-tddft/
zip -r 07-tddft.zip 07-tddft
07-tddft/support.py
+ 84
− 0
View file @ a3340c0b
......@@ -284,6 +284,7 @@ class Volumetric:
k=ijk[:, 2],
opacity=0.5,
hoverinfo="none",
color=color,
)
# THIS IS AN ALTERNATIVE VERSION OF THE GET_MESH METHOD
......@@ -535,11 +536,16 @@ class XYZ:
dynamic_threshold: bool = False,
verbose: bool = False,
) -> None:
if mode == "diffdens":
color = "tomato" if threshold < 0 else "cornflowerblue"
if mode == "mo":
color = "tomato" if threshold < 0 else "cornflowerblue"
self.viewer.add_trace(
Volumetric.from_gbw_file(file, index, mode).get_mesh(
threshold=threshold,
verbose=verbose,
dynamic_threshold=dynamic_threshold,
color=color,
)
)
return
......@@ -807,3 +813,81 @@ class Trajectory:
fig.update_layout(scene_aspectmode="data")
return fig
class Spectrum:
"""
Energies in eV.
"""
units: dict[str, dict[str, str | float | bool]] = {
"ev": {"xaxis_title": "Energy / eV", "factor": 1.0, "inverse": False},
"cm-1": {
"xaxis_title": "Wavenumber / cm^-1",
"factor": 8065.610420,
"inverse": False,
},
"nm": {"xaxis_title": "Wavelength / nm", "factor": 1239.84193, "inverse": True},
}
def __init__(
self, indices: list[int], energies: list[float], intensities: list[float]
):
self.indices = indices
self.energies = energies
self.intensities = intensities
return
@classmethod
def from_tddft_output(cls, output_file: str):
indices = []
energies = []
foscs = []
with open(output_file) as f:
for line in f:
if "ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS" in line:
next(f)
next(f)
next(f)
next(f)
while True:
line = next(f)
if len(line.split()) == 0:
break
line = line.split()
fosc = float(line[3])
if fosc < 0.01:
continue
indices.append(int(line[0]))
energies.append(float(line[1]) / 8065.610420)
foscs.append(fosc)
return cls(indices, energies, foscs)
@staticmethod
def get_impulse_line(
x: list[float], y: list[float]
) -> tuple[list[float | None], list[float | None]]:
new_x: list[float | None] = []
new_y: list[float | None] = []
for xx, yy in zip(x, y):
new_x.append(xx)
new_y.append(0.0)
new_x.append(xx)
new_y.append(yy)
new_x.append(None)
new_y.append(None)
return new_x, new_y
def get_energies_as(self, unit: str) -> list[float]:
unit = unit.lower()
assert unit in ["ev", "cm-1", "nm"], "Unit must be 'ev', 'cm-1' or 'nm'"
if self.units[unit]["inverse"]:
x = [self.units[unit]["factor"] / e for e in self.energies]
else:
x = [self.units[unit]["factor"] * e for e in self.energies]
return x
def get_line_spectrum(
self, unit: str
) -> tuple[list[float | None], list[float | None]]:
return Spectrum.get_impulse_line(self.get_energies_as(unit), self.intensities)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment