Merge branch 'week2' into 'main'

Week2

See merge request !2

02-scf/02-scf.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# B.CTC: SCF\n",
+    "\n",
+    "In diesem Notebook werden wir einen auf das Nötigste reduzierten SCF-Code verwenden, um einen Einblick in die einzelnen Iterationen zu gewinnen."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scf\n",
+    "import numpy as np\n",
+    "from pyscf import gto\n",
+    "import plotly.graph_objects as go"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# To build such a molecule object, specific information is required, which is defined in the following code line. (In this oversimplified case, only the molecular coordinates and the basis set are used.) With this information, PySCF can now calculate the individual system-dependent parameters.\n",
+    "\n",
+    "coordinates='''\n",
+    "    C  0.00  0.00  0.00\n",
+    "    O  0.00  1.54  0.00\n",
+    "    O  0.00 -1.54  0.00\n",
+    "'''\n",
+    "\n",
+    "mol = gto.M(atom=coordinates, basis='sto-3g')\n",
+    "\n",
+    "# nuclei-nuclei potential\n",
+    "E_nuc  = mol.energy_nuc()\n",
+    "\n",
+    "# number of atomic orbitals\n",
+    "NAO = mol.nao\n",
+    "# number of electrons\n",
+    "NEL = mol.nelectron\n",
+    "\n",
+    "# overlap integrals\n",
+    "S = mol.intor('int1e_ovlp')\n",
+    "# kinetic energy integrals\n",
+    "T_e = mol.intor('int1e_kin')\n",
+    "# electrons-nuclei potential energy integrals\n",
+    "V_ne = mol.intor('int1e_nuc')\n",
+    "\n",
+    "# two electron Coulomb integrals\n",
+    "J = mol.intor('int2e')\n",
+    "\n",
+    "\n",
+    "Hc = scf.build_hc(T_e, V_ne)\n",
+    "X = scf.build_ortho_mat(S)\n",
+    "finished_scf = scf.scf(Hc,X,J, E_nuc, NAO, NEL)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = go.Figure()\n",
+    "\n",
+    "fig.add_trace(go.Scatter(\n",
+    "    x = list(range(finished_scf[\"iterations\"])),\n",
+    "    y = finished_scf[\"tot_energies\"],\n",
+    "))\n",
+    "\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = go.Figure()\n",
+    "\n",
+    "orb_eigenvalues = finished_scf[\"orb_eigenvalues\"].T\n",
+    "for i in range(len(orb_eigenvalues)):\n",
+    "    fig.add_trace(go.Scatter(\n",
+    "        x = list(range(finished_scf[\"iterations\"])),\n",
+    "        y = orb_eigenvalues[i],\n",
+    "    ))\n",
+    "\n",
+    "fig.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.set_printoptions(suppress=True, precision=2)\n",
+    "[print((\"{:6.2f}\"*len(line)).format(*line)) for line in finished_scf[\"coefficients\"][-1].T]\n",
+    "# finished_scf[\"coefficients\"][70].T"
+   ]
+  }
+ ],
+ "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
+}

02-scf/scf.py

+#!/usr/bin/env python
+# coding: utf-8
+
+import numpy as np
+
+
+delta_E = 1.0e-8            # Energy convergence criterion
+delta_P = 1.0e-8            # Density matrix convergence criterion
+MAXITER = 100               # Maximum number of SCF Iterations
+
+
+def build_hc(T_e: np.array, V_ne: np.array) -> np.array:
+    return T_e+V_ne
+
+
+def build_ortho_mat(S: np.array) -> np.array:
+    s, U =  np.linalg.eig(S)
+    L    =  np.diag(np.power(s, -0.5))
+    X = np.dot(U,(np.dot(L,np.transpose(U))))
+    return X
+
+
+def build_fock_prime(F: np.array, X: np.array) -> np.array:
+    Fp=np.dot(np.transpose(X),np.dot(F,X))
+    return Fp
+
+
+def eigensort(eigenvalues: np.array, eigenvectors: np.array) -> (np.array, np.array):
+    eigenvectors = eigenvectors[:,eigenvalues.argsort()]
+    eigenvalues = np.sort(eigenvalues)
+    return eigenvalues, eigenvectors
+
+
+def solve_eigenvalue(Fp: np.array, X: np.array) -> (np.array, np.array):
+    e, Cp   =  np.linalg.eigh(Fp)
+    e, Cp   =  eigensort(e, Cp)
+    C       =  np.dot(np.transpose(X),Cp)    
+    return e, C
+
+
+def build_density(C: np.array, NAO: int, NEL: int) -> np.array:
+    P =  np.zeros((NAO, NAO))
+    
+    for i in range(len(C)):
+        for j in range(len(C)):
+            for m in range(int(NEL/2)): #Number of occupied orbitals
+                P[i][j] += C[i][m] * C[j][m]
+    return P
+
+
+def electronic_energy(P: np.array, Hc: np.array, F: np.array, NAO: int) -> float:
+    E = 0
+    for u in range(NAO):
+        for v in range(NAO):
+            E += P[u][v] * (Hc[u][v] + F[u][v])
+    return E
+
+
+def build_new_fock(Hc: np.array, P: np.array, J: np.array, NAO: int) -> np.array:
+    F = np.copy(Hc)
+    for u in range(NAO):
+        for v in range(NAO):
+            for l in range(NAO):
+                for s in range(NAO):
+                    F[u][v] += P[l][s] * ((2 * J[u][v][l][s]) - J[u][l][v][s])
+    return F
+
+
+def compute_rmsd(P0: np.array, P1: np.array) -> float:
+    rmsd = 0
+    for u in range(len(P0)):
+        for v in range(len(P0[u])):
+            rmsd += (P1[u][v] - P0[u][v])**2
+    return np.sqrt(rmsd)
+
+
+def scf(Hc: np.array, X: np.array, J: np.array, E_nuc: float, NAO: int, NEL: int) -> dict:
+    
+    Fp_0 = build_fock_prime(Hc, X)
+    
+    e_0,C_0 = solve_eigenvalue(Fp_0, X) 
+    P_0 = build_density(C_0, NAO, NEL)
+    
+    E_el_0 = electronic_energy(P_0, Hc, Hc, NAO)
+
+    E_tot_0 = E_el_0 + E_nuc
+
+    this_scf = {
+        "tot_energies": [E_tot_0],
+        "orb_eigenvalues": [e_0],
+        "coefficients": [C_0],
+        "iterations": 0
+    }
+
+    # SCF loop
+    for i in range(1,MAXITER+1):
+        
+        F = build_new_fock(Hc, P_0, J, NAO)
+        
+        Fp = build_fock_prime(F, X)
+        
+        e, C = solve_eigenvalue(Fp, X)
+        
+        P = build_density(C, NAO, NEL)
+        
+        E_el = electronic_energy(P, Hc, F, NAO)
+        
+        E_tot = E_el + E_nuc
+
+        this_scf["tot_energies"].append(E_tot)
+        this_scf["orb_eigenvalues"].append(e)
+        this_scf["coefficients"].append(C)
+            
+        if np.abs(E_el-E_el_0)<delta_E:
+            print("E_el converged after: "+str(i)+" Iterations.")
+            print()
+            this_scf["iterations"] = i
+            break
+        
+        E_el_0 = E_el
+                
+        if compute_rmsd(P_0,P)<delta_P:
+            print("P converged after: "+str(i)+" Iterations.")
+            print()
+            this_scf["iterations"] = i
+            break
+        
+        P_0    = np.copy(P)
+    
+    print("   Final RHF results")
+    print("   -----------------")
+    print()
+    print("Total SCF energy = "+str(E_tot))
+    print()
+    print("   Final eigenvalues")
+    print("   -----------------")
+    print()
+    for i in range(len(e)):
+        print(str(i+1)+"  "+str(e[i]))
+
+    this_scf["tot_energies"] = np.array(this_scf["tot_energies"])
+    this_scf["orb_eigenvalues"] = np.array(this_scf["orb_eigenvalues"])
+    
+    return this_scf
+
+
+