CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
howto:newtonx

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
howto:newtonx [2023/08/31 14:05] – [A) Initial conditions and photoabsorption spectra] ahehnhowto:newtonx [2024/01/03 13:09] (current) – oschuett
Line 1: Line 1:
-====== How to run NAMD computations using the CP2K-NEWTONX interface ====== +This page has been moved to: https://manual.cp2k.org/trunk/methods/sampling/newton-x.html
- +
-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 (\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, +
- +
-\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 (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, +
- +
-\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}}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 ] \, . +
-\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 forces +
-  * ''TDDFPT/PRINT/NAMD_PRINT'' with keyword option ''PRINT_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 section 2 of 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''. 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: +
-<code - cp2k.inp> +
-&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 +
-</code> +
-<code - cp2k_vib.inp> +
-&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 +
-</code> +
- +
-The input file ''cp2k.par'' includes all specifications regarding the executable and parallelization setup. +
-<code - cp2k.par> +
- parallel = 16 +
- exec = cp2k.psmp +
-</code> +
- +
-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''. +
- +
-<code - cp2k.par> +
-&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 +
-/ +
-</code> +
- +
-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. +
- +
-<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 +
-</code> +
- +
-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: +
- +
-<code cp2k> +
- 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 +
-</code> +
- +
-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 ===== +
-                                                                             +
howto/newtonx.1693490726.txt.gz · Last modified: by ahehn