{"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.7.8"},"nteract":{"version":"nteract-on-jupyter@2.1.3"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"
\n \n
\n\n# TDSCF with Psi4\n\n*Roberto Di Remigio* ([@robertodr](https://github.com/robertodr))\n\n[European National Competence Centre Sweden](https://enccs.se/)\n \n### **Find this notebook at [https://tinyurl.com/tddft-psicon2020](https://tinyurl.com/tddft-psicon2020)**","metadata":{}},{"cell_type":"markdown","source":"Psi4 1.4 enables you to run TDSCF calculations:\n1. using the random-phase appoximation (RPA). *Credits: Andrew James ([@amjames](https://github.com/amjames)) and Daniel Smith ([@dgasmith](https://https://github.com/dgasmith))*.\n2. using the Tamm-Dancoff approximation (TDA). *Credits: Ruhee D'Cunha ([@rdcunha](https://github.com/rdcunha))*.\n3. from Psithon input files. *Credits: Jeff Schriber ([@jeffschriber](https://github.com/jeffschriber)) and Lori Burns ([@loriab](https://github.com/loriab))*.\n4. writing Python scripts and PsiAPI. *Credits: Andrew James ([@amjames](https://github.com/amjames)) and yours truly*.\n5. including solvent effects with either the polarizable continuum model (PCM) or polarizable embedding (PE). *Credits: Holger Kruse ([@hokru](https://github.com/hokru)) and Maximilian Scheurer ([@maxscheurer](https://github.com/maxscheurer))*.\n6. and easily plot the corresponding spectra. *Credits: yours truly*.\n\n
\n \n \n
\n To run the selected code cell, hit 
Shift + Enter
\n
\n
","metadata":{}},{"cell_type":"markdown","source":"## Using Psithon\n\nThe bridge to Psithon was built by Jeff Schriber ([@jeffschriber](https://github.com/jeffschriber)) and Lori Burns ([@loriab](https://github.com/loriab)).\nA TDSCF calculation can be invoked with two styles:\n- From the `energy` driver\n\n```python\nenergy(\"td-pbe0/cc-pvdz\")\n```\n\n- With a wavefunction object\n\n```python\ne, wfn = energy(\"pbe0/cc-pvdz\", return_wfn=True)\ntdscf(wfn)\n```\n\nWe will ask for 10 roots (states) and to use the Tamm-Dancoff approximation (TDA):\n```python\nset {\n tdscf_states 10\n tdscf_tda true\n}\n```\nother options are detailed in the [manual](http://psicode.org/psi4manual/master/tdscf.html#psithon-keywordshttp://psicode.org/psi4manual/master/tdscf.html#psithon-keywords). The number of roots can be specified either as an integer or as a list. The latter form is mostly useful when dealing with symmetry: you can ask to solve for roots belonging to a specific irrep. This is the complete input file:\n```text\nmolecule moxy {\n0 1\nC 0.152133 -0.035800 0.485797\nC -1.039475 0.615938 -0.061249\nC 1.507144 0.097806 -0.148460\nO -0.828215 -0.788248 -0.239431\nH 0.153725 -0.249258 1.552136\nH -1.863178 0.881921 0.593333\nH -0.949807 1.214210 -0.962771\nH 2.076806 -0.826189 -0.036671\nH 2.074465 0.901788 0.325106\nH 1.414895 0.315852 -1.212218\n}\n\nset_num_threads(4)\n\nset {\n tdscf_states 10\n tdscf_tda true\n}\n\nenergy('td-pbe0/cc-pvdz')\n```\n\nNow we can run with:\n```bash\n$ psi4 -i moxy.inp -o moxy.dat\n```\n\nThe result sections reports excitation energies and transition properties:\n```text\n******************************************************************************************\n********** WARNING **********\n********** Length-gauge rotatory strengths are **NOT** gauge-origin invariant **********\n******************************************************************************************\n\n Excitation Energy Total Energy Oscillator Strength Rotatory Strength \n # Sym: GS->ES (Trans) au eV au au (length) au (velocity) au (length) au (velocity) \n ---- -------------------- --------------- --------------- --------------- --------------- --------------- --------------- ---------------\n 1 A->A (1 A) 0.30939 8.41882 -192.58499 0.0266 0.0185 0.0824 0.0552 \n 2 A->A (1 A) 0.31237 8.50013 -192.58200 0.0052 0.0005 -0.0058 -0.0025 \n 3 A->A (1 A) 0.33179 9.02838 -192.56259 0.0171 0.0077 -0.0356 -0.0294 \n 4 A->A (1 A) 0.33881 9.21956 -192.55556 0.0474 0.0582 -0.0230 -0.0325 \n 5 A->A (1 A) 0.34136 9.28892 -192.55301 0.0169 0.0179 0.0108 0.0181 \n 6 A->A (1 A) 0.34444 9.37258 -192.54994 0.0130 0.0059 -0.0520 -0.0433 \n 7 A->A (1 A) 0.36218 9.85546 -192.53219 0.0231 0.0170 -0.0006 0.0053 \n 8 A->A (1 A) 0.37258 10.13855 -192.52179 0.0307 0.0252 0.0725 0.0708 \n 9 A->A (1 A) 0.37828 10.29347 -192.51610 0.0041 0.0039 0.0001 0.0040 \n 10 A->A (1 A) 0.37904 10.31409 -192.51534 0.1130 0.0521 0.0460 0.0310 \n```\n\nFurther information on the set up of the calculation is reported a bit earlier in the output, together with the progress of the iterative subspace solver:\n```text\n ---------------------------------------------------------\n TDSCF excitation energies \n by Andrew M. James and Daniel G. A. Smith \n ---------------------------------------------------------\n\n ==> Options <==\n\n Residual threshold : 1.0000e-04\n Initial guess : denominators\n Reference : RHF\n Solver type : TDA (Davidson)\n\n\n ==> Requested Excitations <==\n\n 10 singlet states with A symmetry\n\n\n ==> Seeking the lowest 10 singlet states with A symmetry\n\n Generalized Davidson Solver \n By Ruhee Dcunha \n\n ==> Options <==\n\n Max number of iterations = 60 \n Eigenvector tolerance = 1.0000e-04\n Max number of expansion vectors = 2000 \n\n => Iterations <=\n Max[D[value]] Max[|R|] # vectors\n DavidsonSolver iter 1: 3.84154e-01 5.04370e-02 40 \n DavidsonSolver iter 2: 2.40990e-03 3.49018e-02 50 \n DavidsonSolver iter 3: 2.59452e-03 9.67573e-03 60 \n DavidsonSolver iter 4: 1.11697e-04 1.30981e-03 70 \n DavidsonSolver iter 5: 2.39422e-06 1.65681e-04 77 \n DavidsonSolver iter 6: 3.67548e-08 6.20551e-05 78 Converged\n```","metadata":{}},{"cell_type":"markdown","source":"## Using PsiAPI\n\nUsing PsiAPI is as simple as using Psithon. The relevant function is `tdscf_excitations`, which is the same driver called under the hood when running TDSCF from a Psithon input file. In this example, we ran the same example as before, but using the random-phase approximation (RPA). **Note** that the `save_jk` option has to be explicitly set to `True` when working in PsiAPI.","metadata":{}},{"cell_type":"code","source":"import psi4\n\nfrom psi4.driver.procrouting.response.scf_response import tdscf_excitations\n\npsi4.core.set_output_file(\"moxy.out\")\npsi4.core.set_num_threads(4)\n\nmoxy = psi4.geometry(\"\"\"0 1\nC 0.152133 -0.035800 0.485797\nC -1.039475 0.615938 -0.061249\nC 1.507144 0.097806 -0.148460\nO -0.828215 -0.788248 -0.239431\nH 0.153725 -0.249258 1.552136\nH -1.863178 0.881921 0.593333\nH -0.949807 1.214210 -0.962771\nH 2.076806 -0.826189 -0.036671\nH 2.074465 0.901788 0.325106\nH 1.414895 0.315852 -1.212218\n\"\"\", name=\"(S)-methyloxirane\")\n\npsi4.set_options({\n \"save_jk\": True,\n})\n\ne, wfn = psi4.energy(\"HF/cc-pVDZ\", return_wfn=True, molecule=moxy)\nres = tdscf_excitations(wfn, states=10)","metadata":{},"execution_count":1,"outputs":[]},{"cell_type":"markdown","source":"The `res` variable contains all results from the TDSCF calculation. It is a list, with one element per computed root.","metadata":{}},{"cell_type":"code","source":"for k, v in res[0].items():\n print(f\"{k} = {v}\")","metadata":{},"execution_count":2,"outputs":[{"name":"stdout","output_type":"stream","text":"EXCITATION ENERGY = 0.36722797203154167\nELECTRIC DIPOLE TRANSITION MOMENT (LEN) = [-0.00248955 -0.03738262 0.10422942]\nOSCILLATOR STRENGTH (LEN) = 0.0030032956720604273\nELECTRIC DIPOLE TRANSITION MOMENT (VEL) = [ 0.01485039 0.0195916 -0.05823091]\nOSCILLATOR STRENGTH (VEL) = 0.007252903509688956\nMAGNETIC DIPOLE TRANSITION MOMENT = [-0.05400975 0.30059862 0.0395392 ]\nROTATORY STRENGTH (LEN) = -0.006981555278316112\nROTATORY STRENGTH (VEL) = -0.007583130408470832\nSYMMETRY = A\nSPIN = singlet\nRIGHT EIGENVECTOR ALPHA = \nLEFT EIGENVECTOR ALPHA = \nRIGHT EIGENVECTOR BETA = \nLEFT EIGENVECTOR BETA = \n"}]},{"cell_type":"markdown","source":"## Known limitations\n\nThe implementation cannot currently handle the following cases: \n- The use of symmetry with density-fitted two-electron integrals. \n- Excited states of triplet symmetry from a restricted DFT reference. \n- Excited states from an unrestricted reference other than HF or LDA.","metadata":{}},{"cell_type":"markdown","source":"## Theory and implementation\n\nYou can find additional material on response theory in [Norman2011-ad](http://dx.doi.org/10.1039/c1cp21951k) and [Dreuw2005-wp](http://dx.doi.org/10.1021/cr0505627).\n\nThe time-dependent Schrödinger equation is the appropriate\nequation of motion:\n\\begin{equation}\\label{eq:td-schrodinger}\n H|0(t)\\rangle = \\mathrm{i}\\frac{\\partial |0(t)\\rangle}{\\partial t},\n\\end{equation}\nwhere the time-dependent, semiclassical matter-field Hamiltonian is\nused:\n\\begin{equation}\n H = H_0 + V(t).\n\\end{equation}\nThe perturbation is a one-electron operator that is periodic in time $V(t) = V(t+T)$.\nRather than solve for the time-dependent wavefunction, we can target the order-by-order corrections of time-dependent *expectation values*: \n\\begin{equation}\n\\begin{aligned}\n\\langle 0(t) | A | 0(t) \\rangle \n &= \n\\langle 0 | A | 0 \\rangle\n +\n\\sum_{(b, B)} \\langle \\langle A; B \\rangle \\rangle_{\\omega_{b}}\\epsilon_{B}(\\omega_{b})\\exp[-\\mathrm{i}\\omega_{b}t] \\\\\n &+\n \\frac{1}{2}\n\\sum_{(b, B);(c, C)} \\langle \\langle A; B, C \\rangle \\rangle_{\\omega_{b}, \\omega_{c}}\\epsilon_{B}(\\omega_{b})\\epsilon_{C}(\\omega_{c})\\exp[-\\mathrm{i}(\\omega_{b}+\\omega_{c})t]\n + \\ldots\n\\end{aligned}\n\\end{equation}\nThe discrete Fourier coefficients are the **response functions** which we can map to observable properties of the system.\nFor example, the linear response function:\n\\begin{equation}\n\\langle \\langle A; B \\rangle \\rangle_{\\omega}\n=\n\\sum_{i \\neq j} \n\\frac{\\langle i | A | j \\rangle\\langle j | B | i \\rangle}{\\omega - \\omega_{ij}}\n-\n\\frac{\\langle i | B | j \\rangle\\langle j | A | i \\rangle}{\\omega + \\omega_{ij}}\n\\end{equation}\nin the exact eigenbasis of the unperturbed Hamiltonian and it encodes *all* excitation energies of the system (the **poles**) and transition properties (the **residues**).\nThe sum-over-states (SOS) expression requires the whole eigenbasis of the Hamiltonian. However, response properties can be extracted without prior knowledge of a suitable basis, as long as one has access to the molecular electronic Hessian.\n\nIn self-consistent field theories, excitation energies of molecular systems can\nbe obtained as the eigenvalues, $\\omega_{n}$, of the response eigenvalue\nproblem:\n\\begin{equation}\n \\begin{pmatrix}\n \\mathbf{A} & \\mathbf{B} \\\\\n \\mathbf{B}^{*} & \\mathbf{A}^{*}\n \\end{pmatrix}\n \\begin{pmatrix}\n \\mathbf{X}_{n} \\\\\n \\mathbf{Y}_{n} \n \\end{pmatrix}\n =\n \\omega_{n}\n \\begin{pmatrix}\n \\mathbf{1} & \\mathbf{0} \\\\\n \\mathbf{0} & -\\mathbf{1}\n \\end{pmatrix}\n \\begin{pmatrix}\n \\mathbf{X}_{n} \\\\\n \\mathbf{Y}_{n} \n \\end{pmatrix}.\n\\end{equation}\n\nThe $\\mathbf{A}$ and $\\mathbf{B}$ matrices appearing on the left-hand side are\nthe blocks of the molecular electronic\nHessian:\n\\begin{equation}\n\\begin{aligned}\nA_{aibj} &= \\delta_{ij}f_{ab} - \\delta_{ab}f_{ij} + (ai|jb) - \\gamma (ab|ji) + w_{\\mathrm{xc}; aijb} \\\\\nB_{aibj} &= (ai|bj) - \\gamma (aj|ib) + w_{\\mathrm{xc}; aibj} \\\\\n\\end{aligned}\n\\end{equation}\nwhose dimensionality is $(OV)^{2}$, with $O$ and $V$ the number of occupied and\nvirtual molecular orbitals, respectively.\nThis prevents explicit formation of the full Hessian, and subspace iteration\nmethods need to be used to extract the first few roots.\nIn such methods, the eigenvectors are expanded in a subspace of trial vectors,\nwhose dimensionality is greatly lower than that of the full eigenproblem.\nThe Hessian is projected down to this subspace where conventional full\ndiagonalization algorithms can be applied. The subspace is augmented with new\ntrial vectors, until a suitable convergence criterion is met.\nThe efficiency of the subspace solver is determined by the first half-projection\nof the Hessian in the trial subspace, that is, by the efficiency of the routines\nperforming the matrix-vector products.\n\nTThe response eigenvalue equation is **not** an Hermitian\neigenproblem, due to the nonunit metric on the right-hand side. The eigenproblem however has Hamiltonian symmetry: the\nroots appear in pairs $(\\omega_{n}, -\\omega_{n})$, as do the eigenvectors.\nWhen using real orbitals, this symmetry can be preserved by solving an equivalent problem instead. First note that:\n\\begin{equation}\n (\\mathbf{A} - \\mathbf{B})(\\mathbf{A} + \\mathbf{B})| \\mathbf{X}_{n} + \\mathbf{Y}_{n}\\rangle\n =\n \\omega^{2}_{n} | \\mathbf{X}_{n} + \\mathbf{Y}_{n}\\rangle,\n\\end{equation}\nwhich can be brought to the Hermitian form:\n\\begin{equation}\n (\\mathbf{A} - \\mathbf{B})^{\\frac{1}{2}}(\\mathbf{A} + \\mathbf{B})(\\mathbf{A} - \\mathbf{B})^{\\frac{1}{2}} \\mathbf{T}_{n}\n =\n \\omega^{2}_{n} \\mathbf{T}_{n},\n\\end{equation}\nassuming the SCF reference is stable, *i.e.* $(\\mathbf{A}-\\mathbf{B})$ is positive-definite.\nThe paired vectors $| \\mathbf{X}_{n} - \\mathbf{Y}_{n}\\rangle$ are left\neigenvectors and form a biorthonormal set together\nwith the right eigenvectors $| \\mathbf{X}_{n} + \\mathbf{Y}_{n}\\rangle$.\nIn the Tamm-Dancoff approximation (TDA), we solve the Hermitian eigenvalue problem:\n\\begin{equation}\n\\mathbf{A}\\mathbf{X}_{n} = \\omega_{n}\\mathbf{X}_{n},\n\\end{equation}\nusing a regular Davidson solver.","metadata":{}},{"cell_type":"markdown","source":"## Transition properties\n\nThe excitation energies tell only one part of the story. We need to have transition moments to understand whether an excitation is allowed (bright) or not (dark).\n\nThe excitation energies and eigenvectors can be used to compute transition moments, such as electric and magnetic transition dipole moments, and spectroscopic intensities, such as oscillator strengths and rotatory strengths.\nFor example, the transition electric dipole moment in the length gauge is:\n\n\\begin{equation}\n f = \\frac{2}{3} \\omega_{n} \\sum_{u=x,y,z}\\sum_{ia}|(\\mathbf{X}_{n}+\\mathbf{Y}_{n})_{ia}\\mu_{ai, u}|^{2}.\n\\end{equation}\n\nOur implementation calculates the following transition properties:\n- electric dipole moment in the length gauge\n- oscillator strength in the length gauge\n- electric dipole moment in the velocity gauge\n- oscillator strength in the velocity gauge\n- magnetic dipole moment\n- rotatory strength in the length gauge\n\nupon convergence of the solver. It is rather easy to compute new transition properties. The following is a re-implementation of the oscillator strength:","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\n# get dipole moment integrals\nmints = psi4.core.MintsHelper(wfn.basisset())\nC_L = wfn.Ca_subset(\"SO\", \"OCC\")\nC_R = wfn.Ca_subset(\"SO\", \"VIR\")\ndipole = [psi4.core.triplet(C_L, x, C_R, True, False, False) for x in mints.so_dipole()]\n\nfor x in res:\n # Expression in Pedersen, T. B.; Hansen, A. E. Chem. Phys. Lett. 1995, 246, 1\n edtm = np.sqrt(2) * np.array([x[\"RIGHT EIGENVECTOR ALPHA\"].vector_dot(u) for u in dipole])\n f = 2/3 * x[\"EXCITATION ENERGY\"] * np.sum(edtm**2)\n np.testing.assert_allclose(edtm, x[\"ELECTRIC DIPOLE TRANSITION MOMENT (LEN)\"])\n np.testing.assert_allclose(f, x[\"OSCILLATOR STRENGTH (LEN)\"])","metadata":{},"execution_count":3,"outputs":[]},{"cell_type":"markdown","source":"## Plotting spectra\n\nWe can use these observables to plot spectra. We use empirical lineshapes to simulate the broadening observed in experiments.\nFor one-photon absorption (OPA), the function to plot is:\n$$\n\\varepsilon(\\omega) = \n \\frac{4\\pi^{2}N_{\\mathrm{A}}\\omega}{3\\times 1000\\times \\ln(10) (4 \\pi \\epsilon_{0}) n \\hbar c}\n \\sum_{i \\rightarrow j}g_{ij}(\\omega)|\\mathbf{\\mu}_{ij}|^{2}\n$$\nwhile for electronic circular dichroism (ECD):\n$$\n\\Delta\\varepsilon(\\omega) =\n \\frac{16\\pi^{2}N_{\\mathrm{A}}\\omega}{3\\times 1000\\times \\ln(10) (4 \\pi \\epsilon_{0}) n \\hbar c^{2}}\n \\sum_{i \\rightarrow j}g_{ij}(\\omega)\\Im(\\mathbf{\\mu}_{ij}\\cdot\\mathbf{m}_{ij}).\n$$\nSee [Rizzo2011-to](https://doi.org/10.1002/9781118008720.ch2) for details.\n\nI decided to decouple the lineshape fitting from the plotting:\n- We do not entangle Psi4 with any specific plotting library.\n- The user has the freedom to run whatever visualization pipeline.\n\nPsi4 provides the `spectrum` function:\n\n```python\ndef spectrum(*,\n poles: Union[List[float], np.ndarray],\n residues: Union[List[float], np.ndarray],\n kind: str = \"opa\",\n lineshape: str = \"gaussian\",\n gamma: float = 0.2,\n npoints: int = 5000,\n out_units: str = \"nm\") -> Dict[str, np.ndarray]\n```\n\nthat will perform the lineshape fitting of the TDDFT results and return numpy arrays for the $x$ and $y$ values of the plot.\n\nFor both **one-photon absorption** (OPA) and **electronic circular dichroism** (ECD) the *poles* are the excitation energies. The residues are:\n\n- The square of the transition electric dipole moment for OPA: $|\\mu_{ij}|^{2}$.\n- The rotatory strength $R_{ij}$ for ECD. This is the imaginary part of the product of transition electric and magnetic dipole moments: $\\Im(\\mathbf{\\mu}_{ij}\\cdot\\mathbf{m}_{ij})$.\n\nThe implementation assumes these quantities are given in the **length gauge**.\n\nWe can obtain these quantities in atomic units from the list of results returned by the `tdscf_excitations` function.","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\n# get poles and residues to plot OPA and ECD spectra\npoles = [r[\"EXCITATION ENERGY\"] for r in res]\nopa_residues = [np.linalg.norm(r[\"ELECTRIC DIPOLE TRANSITION MOMENT (LEN)\"])**2 for r in res]\necd_residues = [r[\"ROTATORY STRENGTH (LEN)\"] for r in res]","metadata":{},"execution_count":4,"outputs":[]},{"cell_type":"markdown","source":"With these at hand, we invoke the `spectrum` function to perform the convolution with Gaussian lineshapes (inhomogeneous broadening) with empirical linewidth $\\gamma = 0.01\\,\\textrm{a.u.}$ (angular frequency). ","metadata":{}},{"cell_type":"code","source":"from psi4.driver.p4util import spectrum\n\nopa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units=\"nm\")\n\necd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind=\"ECD\", gamma=0.01, out_units=\"nm\")","metadata":{},"execution_count":5,"outputs":[]},{"cell_type":"markdown","source":"All data in input is given in atomic units, but the final spectrum will use nanometers on the $x$ axis and macroscopic units $\\mathrm{L}\\cdot\\mathrm{mol}^{-1}\\cdot\\mathrm{cm}^{-1}$ on the $y$ axis.\n\nThe return value is a dictionary-of-dictionaries with:\n1. The lineshape convolution:\n - `opa_spectrum[\"convolution\"][\"x\"]` a NumPy array with the $x$ axis points of the lineshape convolution.\n - `opa_spectrum[\"convolution\"][\"y\"]` a NumPy array with the $y$ axis points of the lineshape convolution.\n2. The infinitely narrow (stick) representation:\n - `opa_spectrum[\"sticks\"][\"poles\"]` a NumPy array with the poles in the chosen unit of measure.\n - `opa_spectrum[\"sticks\"][\"residues\"]` a NumPy array with the residues in macroscopic units.","metadata":{}},{"cell_type":"code","source":"opa_spectrum[\"sticks\"]","metadata":{},"execution_count":6,"outputs":[{"execution_count":6,"output_type":"execute_result","data":{"text/plain":"{'poles': array([124.07375252, 118.03023459, 115.43690877, 112.12611621,\n 109.32211753, 109.08955351, 106.82893788, 104.23903833,\n 103.73414661, 99.94543184]),\n 'residues': array([3.10116093e+00, 1.07161409e+01, 2.08483761e+02, 1.02133229e+04,\n 5.44012797e+03, 5.80576170e+03, 6.79281405e+03, 1.28870924e+03,\n 1.42292033e+03, 5.19938683e+03])}"},"metadata":{}}]},{"cell_type":"markdown","source":"We can now plot these spectra with any library we want. In this example, I will use [Altair](https://altair-viz.github.io/) and [pandas](https://pandas.pydata.org/). The raw data is stored in NumPy arrays and any plotting library will work!","metadata":{}},{"cell_type":"code","source":"from altair_spectrum import plot_spectrum\n\nopa_plot = plot_spectrum(opa_spectrum,\n title=\"OPA (Gaussian broadening)\",\n x_title=(\"λ\", \"nm\"))\n\necd_plot = plot_spectrum(ecd_spectrum,\n title=\"ECD (Gaussian broadening)\",\n x_title=(\"λ\", \"nm\"),\n y_title=(\"Δε\", \"L⋅mol⁻¹⋅cm⁻¹\"))\n\nopa_plot & ecd_plot","metadata":{},"execution_count":7,"outputs":[{"execution_count":7,"output_type":"execute_result","data":{"text/html":"\n
\n","text/plain":"alt.VConcatChart(...)"},"metadata":{}}]},{"cell_type":"markdown","source":"## Solvent effects\n\nFinally, it is possible to run these calculations including solvent effects, either with the polarizable continuum model (PCM) or with polarizable embedding (PE).","metadata":{}},{"cell_type":"code","source":"import psi4\n\nfrom psi4.driver.procrouting.response.scf_response import tdscf_excitations\n\npsi4.core.set_output_file(\"h2o2+pcm.out\")\npsi4.core.set_num_threads(4)\n\npcm_string = \"\"\"\n Units = Angstrom\n Medium {\n SolverType = IEFPCM\n Solvent = Water\n Nonequilibrium = True\n }\n Cavity {\n Type = GePol\n Area = 1.0\n }\n\"\"\"\n\nh2o2 = psi4.geometry(\"\"\"0 1\nO 0.000000 0.695000 -0.092486 \nO -0.000000 -0.695000 -0.092486 \nH -0.388142 0.895249 0.739888 \nH 0.388142 -0.895249 0.739888 \nsymmetry c1\n\"\"\", name=\"H2O2\")\n\npsi4.set_options({\n \"save_jk\": True,\n \"pcm\": True,\n \"pcm__input\": pcm_string\n})\n\ne, wfn = psi4.energy(\"HF/cc-pVDZ\", return_wfn=True, molecule=h2o2)\nres = tdscf_excitations(wfn, states=10)","metadata":{},"execution_count":8,"outputs":[]},{"cell_type":"markdown","source":"And we can plot the results, of course!","metadata":{}},{"cell_type":"code","source":"import numpy as np\n\n# get poles and residues to plot OPA and ECD spectra\npoles = [r[\"EXCITATION ENERGY\"] for r in res]\nopa_residues = [np.linalg.norm(r[\"ELECTRIC DIPOLE TRANSITION MOMENT (LEN)\"])**2 for r in res]\necd_residues = [r[\"ROTATORY STRENGTH (LEN)\"] for r in res]\n\nfrom psi4.driver.p4util import spectrum\n\nopa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units=\"nm\")\n\necd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind=\"ECD\", gamma=0.01, out_units=\"nm\")\n\nfrom altair_spectrum import plot_spectrum\n\nopa_plot = plot_spectrum(opa_spectrum,\n title=\"OPA (Gaussian broadening)\",\n x_title=(\"λ\", \"nm\"))\n\necd_plot = plot_spectrum(ecd_spectrum,\n title=\"ECD (Gaussian broadening)\",\n x_title=(\"λ\", \"nm\"),\n y_title=(\"Δε\", \"L⋅mol⁻¹⋅cm⁻¹\"))\n\nopa_plot & ecd_plot","metadata":{},"execution_count":9,"outputs":[{"execution_count":9,"output_type":"execute_result","data":{"text/html":"\n
\n","text/plain":"alt.VConcatChart(...)"},"metadata":{}}]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}