From 91d9325da1f83af4bfb71fcbe6af63530c951d51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ole=20Sch=C3=BCtt?= Date: Sun, 10 Dec 2023 22:18:44 +0100 Subject: [PATCH] Manual: Migrate howto:newtonx from wiki (Orig. Author: Anna Hehn) --- docs/methods/sampling/index.md | 2 +- docs/methods/sampling/newton-x.md | 394 ++++++++++++++++++++++++++++++ 2 files changed, 395 insertions(+), 1 deletion(-) create mode 100644 docs/methods/sampling/newton-x.md diff --git a/docs/methods/sampling/index.md b/docs/methods/sampling/index.md index c088120374..ea5ee76bbc 100644 --- a/docs/methods/sampling/index.md +++ b/docs/methods/sampling/index.md @@ -8,7 +8,7 @@ maxdepth: 1 molecular_dynamics metadynamics path_integrals -Surface Hopping with NEWTON-X +newton-x i-PI Gibbs Ensemble Monte Carlo Langevin Dynamics diff --git a/docs/methods/sampling/newton-x.md b/docs/methods/sampling/newton-x.md new file mode 100644 index 0000000000..af83a24124 --- /dev/null +++ b/docs/methods/sampling/newton-x.md @@ -0,0 +1,394 @@ +# 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/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/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