07-tddft.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# B.CTCT: Time-Dependent Density Functional Theory\n",
    "\n",
    "Angeregte Zustände, Näherung, Intensität und Shift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objects as go\n",
    "import subprocess as sub\n",
    "import support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vol = support.Volumetric.from_gbw_file(\"benzene-tddft.gbw\", 3, \"mo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 6*10**-3\n",
    "xyz = support.XYZ.from_xyz_file(\"benzene-opt.xyz\")\n",
    "xyz.get_molecular_viewer()\n",
    "xyz.get_mesh_from_gbw(\"benzene-tddft.gbw\", 0, mode=\"mo\", threshold=threshold)\n",
    "xyz.get_mesh_from_gbw(\"benzene-tddft.gbw\", 0, mode=\"mo\", threshold=-threshold)\n",
    "# xyz.get_mesh_from_cube(\"benzene-tddft.mo5a.cube\", threshold)\n",
    "# xyz.get_mesh_from_cube(\"benzene-tddft.mo5a.cube\", -threshold)\n",
    "xyz.viewer.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tddft_spectrum(output_file: str) -> tuple[list[int], list[float], list[float]]:\n",
    "    indices = []\n",
    "    wavenumbers = []\n",
    "    wavelengths = []\n",
    "    foscs = []\n",
    "    with open(output_file) as f:\n",
    "        for line in f:\n",
    "            if \"ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS\" in line:\n",
    "                next(f)\n",
    "                next(f)\n",
    "                next(f)\n",
    "                next(f)\n",
    "                while True:\n",
    "                    line = next(f)\n",
    "                    if len(line.split()) == 0: break\n",
    "                    line = line.split()\n",
    "                    indices.append(int(line[0]))\n",
    "                    wavenumbers.append(float(line[1]))\n",
    "                    wavelengths.append(float(line[2]))\n",
    "                    foscs.append(float(line[3]))\n",
    "    return indices, wavenumbers, wavelengths, foscs\n",
    "\n",
    "indices, wavenumbers, wavelengths, foscs = get_tddft_spectrum(\"benzene-tddft.out\")\n",
    "fig = go.Figure(layout=go.Layout(\n",
    "    title=\"Benzene TDDFT Spectrum\",\n",
    "    # xaxis_title=\"Wavenumber / cm<sup>-1</sup>\",\n",
    "    xaxis_title=\"Wavelength / nm\",\n",
    "    yaxis_title=\"Oscillator Strength\"))\n",
    "# fig.add_trace(go.Bar(x=wavelengths, y=foscs, text=indices))\n",
    "fig.add_trace(go.Scatter(x=wavelengths, y=foscs, text=indices, mode='text'))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from typing import Iterable\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "class Volumetric:\n",
    "    def __init__(self,\n",
    "                 origin: Iterable[float],\n",
    "                 res_vec1: int, res_vec2: int, res_vec3: int,\n",
    "                 vec1: Iterable[float],\n",
    "                 vec2: Iterable[float],\n",
    "                 vec3: Iterable[float],\n",
    "                 volumetrics: Iterable[float]) -> None:\n",
    "        self.origin = np.array(origin)\n",
    "        self.res_vec1 = res_vec1\n",
    "        self.res_vec2 = res_vec2\n",
    "        self.res_vec3 = res_vec3\n",
    "        self.vec1 = np.array(vec1)\n",
    "        self.vec2 = np.array(vec2)\n",
    "        self.vec3 = np.array(vec3)\n",
    "        self.volumetrics = volumetrics\n",
    "        return\n",
    "\n",
    "    @classmethod\n",
    "    def read_cube_file(cls, file: str):\n",
    "        with open(file, \"r\") as f:\n",
    "            next(f)\n",
    "            next(f)\n",
    "            natoms, x0, y0, z0 = f.readline().split()\n",
    "            res_vec1, x1, y1, z1 = f.readline().split()\n",
    "            res_vec2, x2, y2, z2 = f.readline().split()\n",
    "            res_vec3, x3, y3, z3 = f.readline().split()\n",
    "\n",
    "            natoms, x0, y0, z0 = (abs(int(natoms)), float(x0), float(y0), float(z0))\n",
    "            res_vec1, x1, y1, z1 = (int(res_vec1), float(x1), float(y1), float(z1))\n",
    "            res_vec2, x2, y2, z2 = (int(res_vec2), float(x2), float(y2), float(z2))\n",
    "            res_vec3, x3, y3, z3 = (int(res_vec3), float(x3), float(y3), float(z3))\n",
    "\n",
    "            # skip all the atom coordinates\n",
    "            for i in range(natoms):\n",
    "                next(f)\n",
    "            if \".mo\" in file:\n",
    "                next(f)\n",
    "\n",
    "            volumetrics = f.readlines()\n",
    "        volumetrics = [float(i) for line in volumetrics for i in line.split()]\n",
    "        return cls(\n",
    "            (x0, y0, z0),\n",
    "            res_vec1, res_vec2, res_vec3,\n",
    "            (x1, y1, z1),\n",
    "            (x2, y2, z2),\n",
    "            (x3, y3, z3),\n",
    "            volumetrics\n",
    "        )\n",
    "    \n",
    "    def voxel_generator(self, lower: float, upper: float):\n",
    "        for x in range(self.res_vec1):\n",
    "            for y in range(self.res_vec2):\n",
    "                for z in range(self.res_vec3):\n",
    "                    voxel =  self.volumetrics[z + y * self.res_vec3 + x * (self.res_vec3 * self.res_vec2)]\n",
    "                    if voxel > lower and voxel < upper:\n",
    "                        coord = x * self.vec1 + y * self.vec2 + z * self.vec3\n",
    "                        yield voxel, *coord\n",
    "\n",
    "\n",
    "lower = 6*10**-5\n",
    "vol = Volumetric.read_cube_file(\"benzene-tddft.mo5a.cube\")\n",
    "voxels = [voxel for voxel in vol.voxel_generator(lower, 1)]\n",
    "voxels = np.array(voxels)\n",
    "print(len(voxels))\n",
    "print(voxels[:10])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.cluster as skc\n",
    "\n",
    "cluster = skc.DBSCAN(eps=0.5).fit(voxels[:,1:4])\n",
    "print(cluster.labels_)\n",
    "print(len(cluster.labels_))\n",
    "_ = [print(l, end=\" \") for l in cluster.labels_]\n",
    "print()\n",
    "print(cluster.labels_.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Mesh3d(x=voxels[cluster.labels_ == i][:,1],\n",
    "                            y=voxels[cluster.labels_ == i][:,2],\n",
    "                            z=voxels[cluster.labels_ == i][:,3],\n",
    "                            opacity=0.5,\n",
    "                            alphahull=1\n",
    "                            ))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "for i in range(9):\n",
    "    cluster_i = voxels[cluster.labels_ == i]\n",
    "    fig.add_trace(go.Scatter3d(x=cluster_i[:,1],\n",
    "                            y=cluster_i[:,2],\n",
    "                            z=cluster_i[:,3],\n",
    "                            mode='markers',\n",
    "                            opacity=0.5\n",
    "                            ))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista as pv\n",
    "\n",
    "xyz = np.ndarray((0, 3))\n",
    "ijk = np.ndarray((0, 3))\n",
    "fig = go.Figure()\n",
    "for i in range(9):\n",
    "    cloud = pv.PolyData(voxels[cluster.labels_ == i][:,1:4])\n",
    "    surf = cloud.delaunay_3d(alpha=0.4, tol=0.001, offset=2.5)\n",
    "    surf = surf.extract_geometry()\n",
    "    print(surf.volume)\n",
    "    print(surf.volume)\n",
    "    surf = surf.triangulate()\n",
    "    print(len(surf.points))\n",
    "    print(surf.is_all_triangles)\n",
    "    surf = surf.decimate(0.9)\n",
    "    print(len(surf.points))\n",
    "    ijk = np.concatenate((ijk, surf.faces.reshape(-1, 4)[:, 1:] + len(xyz)))\n",
    "    xyz = np.concatenate((xyz, surf.points))\n",
    "    # xyz = surf.points\n",
    "    # ijk = surf.faces.reshape(-1, 4)[:, 1:]\n",
    "fig.add_trace(\n",
    "    go.Mesh3d(x=xyz[:,0],\n",
    "                    y=xyz[:,1],\n",
    "                    z=xyz[:,2],\n",
    "                    i=ijk[:,0],\n",
    "                    j=ijk[:,1],\n",
    "                    k=ijk[:,2],\n",
    "                    opacity=0.5,\n",
    "                    # alphahull=0\n",
    "    )\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista as pv\n",
    "\n",
    "cloud = pv.PolyData(voxels[:,1:4])\n",
    "surf = cloud.delaunay_3d(alpha=0.4, tol=0.001, offset=2.5)\n",
    "# surf = cloud.reconstruct_surface(nbr_sz=40)\n",
    "# surf = surf.extract_surface()\n",
    "surf = surf.extract_geometry()\n",
    "print(surf.volume)\n",
    "# surf = surf.smooth_taubin(n_iter=100, edge_angle=5, feature_angle=10)\n",
    "print(surf.volume)\n",
    "surf = surf.triangulate()\n",
    "print(len(surf.points))\n",
    "print(surf.is_all_triangles)\n",
    "surf = surf.decimate(0.9)\n",
    "print(len(surf.points))\n",
    "xyz = surf.points\n",
    "ijk = surf.faces.reshape(-1, 4)[:, 1:]\n",
    "# _ = [print(x) for i,x in enumerate(surf.faces) if i % 4 == 0]\n",
    "surf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "fig.add_trace(go.Mesh3d(x=xyz[:,0],\n",
    "                        y=xyz[:,1],\n",
    "                        z=xyz[:,2],\n",
    "                        i=ijk[:,0],\n",
    "                        j=ijk[:,1],\n",
    "                        k=ijk[:,2],\n",
    "                        opacity=0.5,\n",
    "                        # alphahull=0\n",
    "                       ))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}