This is an old revision of the document!
Table of Contents
How to run NAMD computations using the CP2K-NEWTONX interface
This is a short tutorial on how to use the CP2K-NEWTONX interface to a) generate initial conditions to compute photoabsorption spectra and b) to run non-adiabatic dynamics simulations using orbital derivative couplings. A more comprehensive tutorial on all NEWTONX features, including a documentation of the required specifications for the CP2K interface, can be found on the NEWTONX homepage, https://newtonx.org/documentation-tutorials/.
Brief theory recap
The interface enables to use electronic-structure data from CP2K and combine it with the surface hopping module of NEWTONX. Excitation energies $\Omega^M$ and excited-state eigenvectors $\mathbf{X}^M$ to describe the excited state $M$ are provided by CP2K, relying on the Tamm-Dancoff eigenvalue problem,
\begin{equation} \label{tda_equation} \begin{aligned} \mathbf{A} \mathbf{X}^M &= \Omega^M \mathbf{S} \mathbf{X}^M \, , \\ \sum_{\kappa k} [ F_{\mu \kappa \sigma} \delta_{ik} - F_{ik \sigma} S_{\mu \kappa} ] X^M_{\kappa k \sigma} + \sum_{\lambda} K_{\mu \lambda \sigma} [\mathbf{D}^{{\rm{\tiny{X}}}M}] C_{\lambda i \sigma} &= \sum_{\kappa} \Omega^M S_{\mu \kappa} X^M_{\kappa i \sigma} \, , \end{aligned} \end{equation}
with $\mathbf{S}$ representing the conventional atomic-orbital overlap matrix, $\mathbf{F}$ the Kohn-Sham matrix, $\mathbf{K}$ the kernel comprising – depending on the chosen functional – Coulomb, exchange and exchange-correlation contributions, and $\mathbf{C}$ the molecular orbital coefficients. $\mu, \nu, \dots$ denote atomic orbitals, $i, j, \dots$ occupied molecular orbitals. The corresponding excited-state gradient is obtained setting up a variational Lagrangian and taking the derivative with respect to the nuclear coordinates $\mathbf{R}$ (see also https://www.cp2k.org/howto:tddft).
By performing a TDDFPT computation, excitation energies $\Omega^M$, excited-state eigenvectors $\mathbf{X}^M$ and corresponding excited-state gradients $\nabla \Omega^M (\mathbf{R})$ are provided by CP2K. On the so-defined potential energy surfaces, the nuclei are propagated classically relying on the surface hopping code of NEWTONX,
\begin{equation} \label{newtons_eom} \begin{aligned} \mathbf{R}(t + \Delta t) &= \mathbf{R} (t) + \mathbf{v} (t) \Delta t + \frac{1}{2} \mathbf{a}(t) \Delta t^2 \, ,\\ \mathbf{v} (t + \Delta t) &= \mathbf{v} (t) + \frac{1}{2} (\mathbf{a} (t) + \mathbf{a} (t+ \Delta t) ) \Delta t \, , \\ \mathbf{a} (t) &= - \frac{1}{m} \nabla \Omega^M (\mathbf{R}(t)) \, . \end{aligned} \end{equation}
The coefficients $c^M$ of the total wave function $\Psi$ over all excited states $M$ are obtained implying hopping probabilities $P_{M\rightarrow N}$ of Tully's surface hopping,
\begin{equation}\label{surface_hopping} \begin{aligned} \Psi (\mathbf{R}(t)) &= \sum_{M} c^{M} (t) \Psi^M (\mathbf{R}(t)) \\ i \frac{{\rm{d}} c^M (t)}{{\rm{d}} c^M (t)}{{\rm{d}}t} &= \sum_N c^N (t) ( \delta_{MN} E_N (\mathbf{R}(t)) - i \sigma_{MN} (t)) \, , \\ P_{M \rightarrow N} &= {\rm{max}} [ 0, \frac{-2 \Delta t}{| c^M|^2} {\rm{Re}} (c^M c^{N \ast}) \sigma_{MN} ] \, . \end{aligned} \end{equation}
The therefore required non-adiabatic time derivative couplings $\sigma_{MN}$ can be obtained relying on semi-empirical models (Baeck-An; please cite Barbatti et al., Open Research Europe 1, 49 (2021).) or as numerical time derivative couplings (orbital time derivative (OD); please cite Ryabinkin et al., J. Phys. Chem. Lett. 6, 4200 (2015); Barbatti et al., Molecules 21, 1603 (2021).), with the corresponding molecular orbital overlap matrix $\mathbf{S}^{{\rm{\tiny{t-\Delta t,t}}}}$ being provided by CP2K,
\begin{equation}\label{ot_time_deriverative_couplings} \begin{aligned} \sigma_{MN}^{{\rm{\tiny{OD}}}} &= \sum_{ia} X_{ia}^{M} \frac{\partial }{\partial t} X_{ia}^N + \sum_{iab} X_{ia}^M X_{ib}^N S_{ab}^{{\rm{\tiny{t-\Delta t,t}}}} - \sum_{ija} P_{ij} X_{ia}^M X_{ja}^N S_{ji}^{{\rm{\tiny{t-\Delta t,t}}}} \\ S_{pq}^{{\rm{\tiny{t - \Delta t , t}}}} &= \frac{\langle \phi_i (\mathbf{R}(t- \Delta t )) | \phi_j (\mathbf{R} (t)) \rangle}{\Delta t} \, . \end{aligned} \end{equation} $a,b, \dots$ denote virtual molecular orbitals.
General input setup
The input sections for TDDFPT energy and gradient computations are described in the CP2K tutorial https://www.cp2k.org/howto:tddft. To furthermore provide the required CP2K output, subsequently read in by NEWTONX, the following print statements have to be added to the CP2K input files:
FORCE_EVAL/PRINT/FORCES: prints the excited-state forcesTDDFPT/PRINT/NAMD_PRINTwith keyword optionPRINT_PHASES: prints the excited-state eigenvectors in MO format as well as the corresponding phases.VIBRATIONAL_ANALYSIS/PRINT/NAMD_PRINT: prints normal modes to generate initial conditions
It should furthermore be noted that cartesian coordinates have to be provided in terms of the external file “coord.cp2k” and that the number of atoms has to be specified in the CP2K input file in the SUBSYS section.
A) Initial conditions and photoabsorption spectra
The following tutorial to obtain photoabsorption spectra is based on https://vdv.dcf.mybluehost.me/nx/wp-content/uploads/2020/02/tutorial-2_2.pdf.
For the electronic-structure calculation with CP2K, a cp2k.inp and cp2k.par file as well as a coordinate file named coord.cp2k has to be provided in a subdirectory called JOB_AD, with cp2k.inp including all required print sections stated above. Furthermore, to generate the initial conditions, the initqp_input file requires to specify iprog = 10 for CP2K and file_nmodes = cp2k.eig to refer to the corresponding output file comprising the normal modes provided by CP2K. All other keywords are to be chosen as outlined in the corresponding NEWTONX tutorial.
Examplary input files for computing the absorption spectrum of a water molecule are given below:
- h2o_cp2k.inp
The resulting output file of the initcond.pl script of NEWTONX states that the read-in cartesian normal modes are first transfered to mass-weighted normal modes. <code cp2k> Cartesian normal modes (1/sqrt(amu)) 0.00 0.00 0.00 0.00 0.00 0.00 1523.92 3851.12 0.0000 -0.0492 0.0001 -0.1268 0.5632 -0.0083 0.0000 -0.0000 -0.0886 0.0000 -0.0000 -0.0169 0.0047 0.5777 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.5630 0.1269 0.0155 -0.0715 0.0487 0.0001 0.3905 -0.0004 -0.1267 0.5632 -0.0082 -0.4184 -0.5910 0.7043 0.0008 0.7071 -0.0162 0.0040 0.5768 0.0000 0.0000 -0.0001 -0.5885 0.0007 0.5630 0.1270 0.0155 0.5678 -0.3867 0.0000 0.3905 -0.0004 -0.1267 0.5632 -0.0083 0.4184 0.5910 0.7043 -0.0009 -0.7071 -0.0170 0.0051 0.5768 0.0000 0.0000 -0.0000 0.5885 -0.0007 0.5630 0.1269 0.0154 0.5678 -0.3867 3986.44 0.0712 -0.0000 0.0000 -0.5650 0.0000 -0.4222 -0.5650 0.0000 0.4222 Mass weighted normal modes Frequencies will be multiplied by ANH_F = 1.00000 0.00 0.00 0.00 0.00 0.00 0.00 1523.92 3851.12 0.0001 -0.1967 0.0006 -0.5069 2.2526 -0.0330 0.0000 -0.0000 -0.3543 0.0000 -0.0000 -0.0677 0.0186 2.3104 0.0000 -0.0000 -0.0001 -0.0000 -0.0002 2.2517 0.5077 0.0619 -0.2861 0.1949 0.0001 0.3920 -0.0004 -0.1272 0.5654 -0.0083 -0.4200 -0.5933 0.7071 0.0008 0.7099 -0.0162 0.0040 0.5791 0.0000 0.0000 -0.0001 -0.5908 0.0007 0.5652 0.1275 0.0155 0.5700 -0.3882 0.0000 0.3921 -0.0004 -0.1272 0.5654 -0.0083 0.4200 0.5933 0.7071 -0.0009 -0.7099 -0.0171 0.0051 0.5790 0.0000 0.0000 -0.0000 0.5908 -0.0007 0.5652 0.1274 0.0155 0.5700 -0.3882 3986.44 0.2847 -0.0000 0.0000 -0.5672 0.0000 -0.4238 -0.5672 0.0000 0.4238
The thereon based initial conditions are summarized in the output files dubbed “final_output”, comprising geometries and velocities, as examplarily given below,
Initial condition = 1
Geometry in COLUMBUS and NX input format:
o 8.0 5.00630777 5.00000001 4.46399957 15.99491464
h 1.0 6.37684065 5.00000128 5.50815661 1.00782504
h 1.0 3.52303474 5.00000149 5.58297278 1.00782504
Velocity in NX input format:
-0.000089112 0.000000000 -0.000020915
0.000417197 0.000000002 0.000694479
0.000997296 0.000000013 -0.000362483
Epot of initial state (eV): 0.0865 Epot of final state (eV): 19.0799
Vertical excitation (eV): 18.9935 Is Ev in the required range? YES
Ekin of initial state (eV): 0.0479 Etot of initial state (eV): 0.1343
Oscillator strength: 0.1221
State: 10
Moreover, the output file cross-section.dat comprises the data points of the computing photoabsorption spectrum, as shown below.