# Surface Hopping with NEWTON-X 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, . ## 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, $$ \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} \, , $$ 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 [](../properties/optical/tddft)). By performing a TDDFPT computation, excitation energies $\Omega^M (\mathbf{R}(t))$, excited-state eigenvectors $\mathbf{X}^M (\mathbf{R}(t))$ and corresponding excited-state gradients $\nabla \Omega^M (\mathbf{R}(t))$ are provided by CP2K. On the so-defined potential energy surfaces, the nuclei are propagated classically relying on the surface hopping code of NEWTONX, $$ \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)) \, . $$ The coefficients $c^M (t)$ of the total wave function $\Psi (\mathbf{R}(t))$ over all excited states $M$ are obtained implying hopping probabilities $P_{M\rightarrow N}$ of Tully's surface hopping, $$ \Psi (\mathbf{R}(t)) &= \sum_{M} c^{M} (t) \Psi^M (\mathbf{R}(t)) \\ i \frac{{\rm{d}} c^M (t)}{{\rm{d}}t} &= \sum_N c^N (t) \left ( \delta_{MN} E_N (\mathbf{R}(t)) - i \sigma_{MN} (t) \right ) \, , \\ P_{M \rightarrow N} &= {\rm{max}} \left [ 0, \frac{-2 \Delta t}{| c^M|^2} {\rm{Re}} (c^M c^{N \ast}) \sigma_{MN} \right ] \, . $$ 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)](https://doi.org/10.12688/openreseurope.13624.1).) or as numerical time derivative couplings (orbital time derivative (OD); please cite [Ryabinkin et al., J. Phys. Chem. Lett. 6, 4200 (2015)](https://doi.org/10.1021/acs.jpclett.5b02062); [Barbatti et al., Molecules 21, 1603 (2016)](https://doi.org/10.3390/molecules21111603).), with the corresponding molecular orbital overlap matrix $\mathbf{S}^{{\rm{\tiny{t-\Delta t,t}}}}$ being provided by CP2K, $$ \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} \, . $$ $a,b, \dots$ denote virtual molecular orbitals. ## General input setup The input sections for TDDFPT energy and gradient computations are described in [](../properties/optical/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](#CP2K_INPUT.FORCE_EVAL.PRINT.FORCES): prints the excited-state forces - [TDDFPT.PRINT.NAMD_PRINT](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.TDDFPT.PRINT.NAMD_PRINT) with keyword option [PRINT_PHASES](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.TDDFPT.PRINT.NAMD_PRINT.PRINT_PHASES): prints the excited-state eigenvectors in MO format as well as the corresponding phases. - [VIBRATIONAL_ANALYSIS.PRINT.NAMD_PRINT](#CP2K_INPUT.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](#CP2K_INPUT.FORCE_EVAL.SUBSYS) section. ## A) Initial conditions and photoabsorption spectra The following tutorial to obtain photoabsorption spectra is based on section 2 of . 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`. Furthermore, a vibrational analysis computation has to be performed to provide cartesian normal modes, with the input file including the corresponding `NAMD print` section. Examplary input files for computing the absorption spectrum as well as for performing a vibrational analysis for a single water molecule with CP2K are given below: ```none &GLOBAL PROJECT excited_states_for_h2o RUN_TYPE ENERGY PREFERRED_DIAG_LIBRARY SL PRINT_LEVEL medium &END GLOBAL &FORCE_EVAL &PRINT # print statement for ground-state or excited-state forces &FORCES &END FORCES &END PRINT METHOD Quickstep &PROPERTIES &TDDFPT # TDDFPT input section to compute 10 excited states &DIPOLE_MOMENTS DIPOLE_FORM LENGTH &END DIPOLE_MOMENTS KERNEL FULL NSTATES 10 MAX_ITER 100 MAX_KV 20 CONVERGENCE [eV] 1.0e-5 RKS_TRIPLETS F &PRINT # NAMD print section to print excited-state eigenvectors &NAMD_PRINT PRINT_VIRTUALS T PRINT_PHASES T &END NAMD_PRINT &END PRINT &END TDDFPT &END PROPERTIES &DFT &QS METHOD GAPW EPS_DEFAULT 1.0E-17 &END QS &SCF SCF_GUESS restart &OT PRECONDITIONER FULL_ALL MINIMIZER DIIS &END OT &OUTER_SCF MAX_SCF 900 EPS_SCF 1.0E-7 &END OUTER_SCF MAX_SCF 10 EPS_SCF 1.0E-7 &END SCF POTENTIAL_FILE_NAME POTENTIAL BASIS_SET_FILE_NAME EMSL_BASIS_SETS &MGRID CUTOFF 1000 REL_CUTOFF 100 NGRIDS 5 &END MGRID &POISSON PERIODIC NONE PSOLVER MT &END &XC &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC &END DFT &SUBSYS &CELL ABC 8.0 8.0 8.0 PERIODIC NONE &END CELL # Coordinates are provided externally for the interface &COORD @include coord.cp2k &END COORD &TOPOLOGY &CENTER_COORDINATES T &END NATOMS 3 # specifying number of atoms for NEWTONX CONNECTIVITY OFF &END TOPOLOGY &KIND H BASIS_SET 6-311Gxx POTENTIAL ALL &END KIND &KIND O BASIS_SET 6-311Gxx POTENTIAL ALL &END KIND &END SUBSYS &END FORCE_EVAL ``` ```none &GLOBAL PROJECT normal_modes_for_h2o RUN_TYPE VIBRATIONAL_ANALYSIS #computing normal modes to generate initial conditions PREFERRED_DIAG_LIBRARY SL PRINT_LEVEL medium &END GLOBAL &FORCE_EVAL &PRINT &FORCES &END FORCES &END PRINT METHOD Quickstep &DFT &QS METHOD GAPW # GAPW enables comparison with all-electron molecular program codes like Turbomole EPS_DEFAULT 1.0E-17 &END QS &SCF SCF_GUESS restart &OT PRECONDITIONER FULL_ALL MINIMIZER DIIS &END OT &OUTER_SCF MAX_SCF 900 EPS_SCF 1.0E-7 &END OUTER_SCF MAX_SCF 10 EPS_SCF 1.0E-7 &END SCF POTENTIAL_FILE_NAME POTENTIAL BASIS_SET_FILE_NAME EMSL_BASIS_SETS &MGRID CUTOFF 1000 REL_CUTOFF 100 NGRIDS 5 &END MGRID &POISSON PERIODIC NONE PSOLVER MT &END &XC &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC &END DFT &SUBSYS &CELL ABC 8.0 8.0 8.0 PERIODIC NONE &END CELL # coordinates must be provided as external file for NEWTONX &COORD @include coord.cp2k &END COORD &TOPOLOGY &CENTER_COORDINATES T &END NATOMS 3 CONNECTIVITY OFF &END TOPOLOGY &KIND H BASIS_SET 6-311Gxx POTENTIAL ALL &END KIND &KIND O BASIS_SET 6-311Gxx POTENTIAL ALL &END KIND &END SUBSYS &END FORCE_EVAL &VIBRATIONAL_ANALYSIS &PRINT &NAMD_PRINT # keyword to enable printing of cartesian normal modes &END NAMD_PRINT &END PRINT DX 0.001 &END VIBRATIONAL_ANALYSIS ``` The input file `cp2k.par` includes all specifications regarding the executable and parallelization setup. ```none parallel = 16 exec = cp2k.psmp ``` Furthermore, a `initqp_input` file has to be generated for NEWTONX following the instructions given in the NEWTONX tutorial. Specifications for CP2K in the `initqp_input` file are the following: - The file comprising the normal modes of the CP2K frequency computation -- for the above input provided as `normal_modes_for_h2o-VIBRATIONS-1.eig`-- has to be specified as `file_nmodes = normal_modes_for_h2o-VIBRATIONS-1.eig`. - The electronic structure program has to be specified as CP2K by defining `iprog = 10`. ```none &dat nact = 2 iprog = 10 numat = 3 npoints = 500 file_geom = geom file_nmodes = normal_modes_for_h2o-VIBRATIONS-1.eig anh_f = 1 rescale = n temp = 0 ics_flg = n chk_e = 1 nis = 1 nfs = 11 kvert = 1 de = 100 prog = 14 iseed = 0 lvprt = 1 / ``` After providing the excited-state CP2K computation based on input file `h2o_cp2k.inp` in the subdirectory `JOB_AD`, the normal modes `normal_modes_for_h2o-VIBRATIONS-1.eig` of the frequency computation and the `initqp_input` file for NEWTONX, the script initcond.pl of NEWTONX can be executed to generate initial conditions. The resulting initcond-output file of NEWTONX, it is first stated that the read-in cartesian normal modes are transferred to mass-weighted normal modes. ```none 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 external output files for each state, dubbed "final_output_XXX", comprising information on the various geometries and velocities as examplarily given below: ```none 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 ``` Based on the initial conditions, the broadened photoabsorption spectrum can be computed with the nxinp script. As outlined in section 2.7 of the cited NEWTONX tutorial, the so-obtained output file `cross-section.dat` comprises the data points of the computed photoabsorption spectrum as visualized below: ## B) Non-adiabatic dynamics using orbital determinant derivatives