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

Draft on week 7 with mo viewer and stuff

parent 5d6b069f
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
07-tddft/07-tddft.ipynb 0 → 100644
+ 337
− 0
View file @ 2c0c2e87
%% Cell type:markdown id: tags:
# B.CTCT: Time-Dependent Density Functional Theory
Angeregte Zustände, Näherung, Intensität und Shift
%% Cell type:markdown id: tags:
%% 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")
```
%% Cell type:code 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()
```
%% 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
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()
```
%% Cell type:code id: tags:
``` python
from typing import Iterable
import numpy as np
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])
```
%% Cell type:code id: tags:
``` python
import sklearn.cluster as skc
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())
```
%% 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()
```
%% Cell type:code id: tags:
``` 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:
``` python
import pyvista as pv
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()
```
%% 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()
```
%% Cell type:code 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()
```
%% Cell type:code id: tags:
``` python
```
07-tddft/atom_properties.py 0 → 100644
+ 292
− 0
View file @ 2c0c2e87
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 27 11:51:18 2024
@author: benjaminlear
Adapted by M.Buchhorn, thank you Ben!
"""
atom_colors = [
"#000000", # 0 Unknown
"#dddddd", # 1 H
"cyan", # 2 He
"violet", # 3 Li
"#B22222", # 4 Be
"beige", # 5 B
"#444444", # 6 C
"blue", # 7 N
"red", # 8 O
"green", # 9 F
"cyan", # 10 Ne
"#FF1493", # 11 Na
"#A52A2A", # 12 Mg
"#800080", # 13 Al
"#D2691E", # 14 Si
"#228B22", # 15 P
"yellow", # 16 S
"green", # 17 Cl
"cyan", # 18 Ar
"#0000CD", # 19 K
"#8B4513", # 20 Ca
"#8FBC8F", # 21 Sc
"#B8860B", # 22 Ti
"#4682B4", # 23 V
"#B22222", # 24 Cr
"#2E8B57", # 25 Mn
"#FFD700", # 26 Fe
"#DAA520", # 27 Co
"#A52A2A", # 28 Ni
"#4169E1", # 29 Cu
"#708090", # 30 Zn
"#C0C0C0", # 31 Ga
"#808000", # 32 Ge
"#00FF00", # 33 As
"#00CED1", # 34 Se
"#0000FF", # 35 Br
"cyan", # 36 Kr
"#8A2BE2", # 37 Rb
"#8B4513", # 38 Sr
"#9400D3", # 39 Y
"#FF4500", # 40 Zr
"#4682B4", # 41 Nb
"#B22222", # 42 Mo
"#FFD700", # 43 Tc
"#A52A2A", # 44 Ru
"#228B22", # 45 Rh
"#8B4513", # 46 Pd
"#00CED1", # 47 Ag
"#696969", # 48 Cd
"#C0C0C0", # 49 In
"#808000", # 50 Sn
"#228B22", # 51 Sb
"#D2691E", # 52 Te
"#00FF00", # 53 I
"cyan", # 54 Xe
"#8A2BE2", # 55 Cs
"#8B4513", # 56 Ba
"#FF4500", # 57 La
"#FF1493", # 58 Ce
"#9400D3", # 59 Pr
"#4682B4", # 60 Nd
"#DAA520", # 61 Pm
"#B22222", # 62 Sm
"#228B22", # 63 Eu
"#A52A2A", # 64 Gd
"#B8860B", # 65 Tb
"#8FBC8F", # 66 Dy
"#4682B4", # 67 Ho
"#DAA520", # 68 Er
"#B22222", # 69 Tm
"#228B22", # 70 Yb
"#FF4500", # 71 Lu
"#A52A2A", # 72 Hf
"#4682B4", # 73 Ta
"#DAA520", # 74 W
]
atom_symbols = [
"unk", # 0
"H", # 1
"He", # 2
"Li", # 3
"Be", # 4
"B", # 5
"C", # 6
"N", # 7
"O", # 8
"F", # 9
"Ne", # 10
"Na", # 11
"Mg", # 12
"Al", # 13
"Si", # 14
"P", # 15
"S", # 16
"Cl", # 17
"Ar", # 18
"K", # 19
"Ca", # 20
"Sc", # 21
"Ti", # 22
"V", # 23
"Cr", # 24
"Mn", # 25
"Fe", # 26
"Co", # 27
"Ni", # 28
"Cu", # 29
"Zn", # 30
"Ga", # 31
"Ge", # 32
"As", # 33
"Se", # 34
"Br", # 35
"Kr", # 36
"Rb", # 37
"Sr", # 38
"Y", # 39
"Zr", # 40
"Nb", # 41
"Mo", # 42
"Tc", # 43
"Ru", # 44
"Rh", # 45
"Pd", # 46
"Ag", # 47
"Cd", # 48
"In", # 49
"Sn", # 50
"Sb", # 51
"Te", # 52
"I", # 53
"Xe", # 54
"Cs", # 55
"Ba", # 56
"La", # 57
"Ce", # 58
"Pr", # 59
"Nd", # 60
"Pm", # 61
"Sm", # 62
"Eu", # 63
"Gd", # 64
"Tb", # 65
"Dy", # 66
"Ho", # 67
"Er", # 68
"Tm", # 69
"Yb", # 70
"Lu", # 71
"Hf", # 72
"Ta", # 73
"W", # 74
"Re", # 75
"Os", # 76
"Ir", # 77
"Pt", # 78
"Au", # 79
"Hg", # 80
"Tl", # 81
"Pb", # 82
"Bi", # 83
"Po", # 84
"At", # 85
"Rn", # 86
"Fr", # 87
"Ra", # 88
"Ac", # 89
"Th", # 90
"Pa", # 91
"U", # 92
"Np", # 93
"Pu", # 94
"Am", # 95
"Cm", # 96
"Bk", # 97
"Cf", # 98
"Es", # 99
"Fm", # 100
"Md", # 101
"No", # 102
"Lr", # 103
"Rf", # 104
"Db", # 105
"Sg", # 106
"Bh", # 107
"Hs", # 108
"Mt", # 109
"Ds", # 110
"Rg", # 111
"Cn", # 112
"Nh", # 113
"Fl", # 114
"Mc", # 115
"Lv", # 116
"Ts", # 117
"Og", # 118
]
vdw_radii = [
1.70, # 0 Unknown, same as Carbon
1.20, # 1 H
1.40, # 2 He
1.82, # 3 Li
1.53, # 4 Be
1.92, # 5 B
1.70, # 6 C
1.55, # 7 N
1.52, # 8 O
1.47, # 9 F
1.54, # 10 Ne
2.27, # 11 Na
1.73, # 12 Mg
1.84, # 13 Al
2.10, # 14 Si
1.80, # 15 P
1.80, # 16 S
1.75, # 17 Cl
1.88, # 18 Ar
2.75, # 19 K
2.31, # 20 Ca
2.11, # 21 Sc
2.00, # 22 Ti
2.00, # 23 V
2.00, # 24 Cr
2.00, # 25 Mn
2.00, # 26 Fe
2.00, # 27 Co
1.63, # 28 Ni
1.40, # 29 Cu
1.39, # 30 Zn
1.87, # 31 Ga
2.11, # 32 Ge
1.85, # 33 As
1.90, # 34 Se
1.85, # 35 Br
2.02, # 36 Kr
3.03, # 37 Rb
2.49, # 38 Sr
2.19, # 39 Y
2.06, # 40 Zr
2.00, # 41 Nb
2.00, # 42 Mo
2.00, # 43 Tc
2.00, # 44 Ru
2.00, # 45 Rh
2.00, # 46 Pd
1.72, # 47 Ag
1.58, # 48 Cd
1.93, # 49 In
2.17, # 50 Sn
2.00, # 51 Sb
2.06, # 52 Te
1.98, # 53 I
2.16, # 54 Xe
3.43, # 55 Cs
2.68, # 56 Ba
2.50, # 57 La
2.48, # 58 Ce
2.47, # 59 Pr
2.45, # 60 Nd
2.43, # 61 Pm
2.42, # 62 Sm
2.40, # 63 Eu
2.38, # 64 Gd
2.37, # 65 Tb
2.35, # 66 Dy
2.33, # 67 Ho
2.32, # 68 Er
2.30, # 69 Tm
2.28, # 70 Yb
2.27, # 71 Lu
2.16, # 72 Hf
2.09, # 73 Ta
2.02, # 74 W
]
vdw_radii_dict = dict(zip(atom_symbols, vdw_radii))
atom_colors_dict = dict(zip(atom_symbols, atom_colors))
\ No newline at end of file
07-tddft/qc_parser.py 0 → 100644
+ 528
− 0
View file @ 2c0c2e87
This diff is collapsed.
07-tddft/qc_support.py 0 → 100644
+ 118
− 0
View file @ 2c0c2e87
"""
qc_support.py
Support functions for Quantum Chemistry (QC) analysis and visualization.
This module provides utility functions to support QC analysis and visualization,
including 3D molecular visualization and file conversion.
Meant to be used with Jupyter Notebooks in the computational chemistry course
of the TU Darmstadt.
"""
import py3Dmol as pm
from pathlib import Path
import subprocess as sub
import os
from ipywidgets import interact,fixed,IntSlider
import ipywidgets
def mol_to_xyz(file_path: str) -> None:
"""
Convert a mol file to an xyz file. A new file with the same name
but .xyz file ending is created.
@param file_path: name/path of the mol file that is converted.
@return: None
"""
lines = Path(file_path).read_text().split("\n")
natoms = int(lines[3].split()[0])
lines = lines[4:4+natoms]
lines = [line.split()[:4] for line in lines]
lines = [f"{line[3]} {line[0]} {line[1]} {line[2]}" for line in lines]
out = [f"{natoms}", f"converted from {file_path}"]
out += lines
Path(file_path.replace(".mol", ".xyz")).write_text("\n".join(out))
def show_trajectory(file_path: str) -> pm.view:
v = pm.view()
v.addModelsAsFrames(Path(file_path).read_text(), "xyz")
v.zoomTo()
v.setStyle({'sphere':{'radius':0.2},'stick':{'radius':0.1}})
v.animate({"loop" : "forward", "interval" : "2000"})
return v
def show_mol(file_path: str, numbers: bool=False, labels: list[str]=[]) -> pm.view:
v = pm.view()
xyz = Path(file_path).read_text()
v.addModel(xyz, "xyz")
v.zoomTo()
v.setStyle({'sphere':{'radius':0.2},'stick':{'radius':0.1}})
if numbers:
num = int(xyz.split()[0])
for i in range(num):
v.addLabel(f"{i}",
{"backgroundColor" : "white",
"backgroundOpacity" : "0.2",
"fontColor" : "black"},
{"serial" : f"{i}"})
if len(labels) != 0:
num = int(xyz.split()[0])
assert len(labels) == num, "Length of label list does not match number of atoms."
for i,label in enumerate(labels):
v.addLabel(f"{label}",
{"backgroundColor" : "white",
"backgroundOpacity" : "0.2",
"fontColor" : "black"},
{"serial" : f"{i}"})
return v
def read_xyz_trajectory_str(file: str) -> list[str]:
with open(file) as f:
lines = f.readlines()
block_size = int(lines[0]) +2
frames = []
for i in range(0,len(lines), block_size):
block = lines[i:i+block_size]
frames.append("".join(block))
return frames
def show_vibration(file_path: str, index: int) -> pm.view:
"""
Returns a 3Dmol view that is an animation of a vibrational mode.
@param file_path: String. Path to an ORCA output file (.out) or hessian (.hess)
@param index: Integer. Number of the vibrational mode to be animated
@return: 3Dmol view. Animation of the given vibration.
"""
assert ".out" in file_path or ".hess" in file_path, f"""\
The given file is not compatible. Provide either an ORCA output
file (.out) or a hessian (.hess) file.
"""
proc = sub.run(f"orca_pltvib {file_path} {index}", shell=True, text=True, capture_output=True)
v = show_trajectory(f"{file_path}.v{index:03d}.xyz")
os.remove(f"{file_path}.v{index:03d}.xyz")
return v
class FrameViewer():
def __init__(self, traj_file: str) -> None:
self.frames = read_xyz_trajectory_str(traj_file)
self.v = pm.view()
self.v.addModel(self.frames[0], "xyz")
self.v.zoomTo()
self.v.setStyle({'sphere':{'radius':0.2},'stick':{'radius':0.1}})
return
def show_frame(self, index=0):
"""Updates the FrameViewer to show the structure of the
molecule at the given index."""
self.v.removeAllModels()
self.v.addModel(self.frames[index], "xyz")
self.v.setStyle({'sphere':{'radius':0.2},'stick':{'radius':0.1}})
return self.v.show()
def frame_slider(self):
"""Returns a Slider with a view on the given optimization
trajectory. The slider can be used to change the frame."""
return interact(self.show_frame,
index=ipywidgets.IntSlider(min=0,max=len(self.frames)-1, step=1))
\ No newline at end of file
07-tddft/support.py 0 → 100644
+ 637
− 0
View file @ 2c0c2e87
This diff is collapsed.
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