4.2.49. RHODYN

The program RHODYN is used for the density matrix based time-dependent configuration interaction (\(\rho\)-TDCI) calculations. The dynamics of a system initiated by electromagnetic pulse can be obtained by means of density matrix propagation with the RASSI Hamiltonian [109]. Interaction with the field is included in dipole approximation \(H_{\text{int}} = - \vec{d} \cdot \vec{E}\). The theory and methodology can be found in [136]. The feature of the method is the treating core states and valence states on the same footing if the Hamiltonian is chosen properly; depending on the chosen active space different processes can be considered. Because of the faster timescale (atto- to femtosecond), core-states dynamics are even more adequately described in the framework of TDCI since electron dynamics is isolated from nuclear effects. Nonetheless, such effects as population and phase relaxation due to nuclear bath, ionization, and Auger decay can be taken into account (see keywords). Output of the program:

  • Populations of basis states

    Configuration state functions, spin-free, spin–orbit coupled states can be chosen as a basis.

  • Full density matrix

    The full matrix can be output after a given timestep. It can be used to compute various time-dependent observables.

  • Time-evolution of dipole moment

    The value of an observable can be accessed as \(\langle\vec{d}\rangle = \Tr[\rho \vec{d}]\). If requested it is written to file DIPOLE.dat.

These data can be used to compute

  • Linear and non-linear spectra, e.g., TD absorption, emission spectra, high harmonic generation, etc. Dependencies

The program workflow requires some input files to be present with wave function specifications from RASSCF, CASPT2 and their properties from RASSI in HDF5 format. All these dependencies, however, are totally transparent to the user. Files Input files


HDF5 file generated by RASSCF or CASPT2 modules. X stands for number of file. Number of files should be equal to the number of spin manifolds stated with keyword NRSM.


RASSI output file in HDF5 format, calculated on given RASSCF wave functions. It should contain also additional sets (arrays) in HDF5 file which are not standardly included. Writing these sets to HDF5 file is activated with the keyword RHODyn in RASSI.


File with Huang–Rhys factors. If keyword KEXT is active then this file has to contain dissipation rate \(k\) matrix. This file should be generated by the user. In general, it provides a 3-rank tensor, see [137] Eq. 8, which is stored in molden format. Section [FREQ] should contain frequencies of normal modes and section [HR-FACT] gives a matrix (or line) of HR factors for each normal mode.


File with initial density matrix (as an \(N\times N\) ASCII matrix). When keyword POPUlation style set to FROMFILE, then this file should be present. This file should be generated by the user. Output files


Stores the applied electric field. Duplicate in file RDOUT.


This file contains the diagonal density matrix in CSF basis.


This file contains the diagonal density matrix in SF basis.


This file contains the diagonal density matrix in SO basis


Intermediate file in HDF5 format which is automatically created and contains all required ingredients for propagation: transformation matrices from one basis to another, full and spin–orbit coupling Hamiltonians, transition dipole matrix, initial density matrix, and Dyson amplitudes. Once created this file can be used as an input file with keyword RUNMode set to 2.


Main output file. Here, almost all relevant information is stored in HDF5 format. By far, it includes the electric field data, Hamiltonian used for propagation, decay matrix. The time-evolution of the DM diagonal is also stored here; thus it is redundant. Optionally, full density matrix and emission spectra can be written.


This file stores the TD-dipole moment data Input General keywords


The number of spin manifolds \(n\)


The total number of states involved in the dynamics followed by a list of relevant states. If all states from input files are supposed to be used this list can be replaced by the word ALL.


This keyword should be followed by \(n\) lines (math:n is the number of spin manifolds, see NRSM), each containing four numbers: number of determinants, number of CSFs, number of roots, and spin multiplicity for each spin manifold. See examples below.

POPUlation style

Describes how the initial DM is to be constructed depending on the chosen basis. Available options: CSF, SF, SF_THERMAL, SO, SO_THERMAL, FROMFILE.

CSF — means that given CSFs (specified by keyword NRPO) are populated evenly.

SF — same as CSF but for spin-free states.

SF_THERMAL — SFs will be populated according to the Boltzmann distribution for a given temperature, see keyword TEMP.

SO — same as SF and CSF but for spin–orbit-coupled basis states.

SO_THERMAL — same as SF_THERMAL but for SO basis.

FROMFILE — read the full initial DM from file INDENS.


Numbers of states to be populated if POPU is set to CSF, SF, or SO. The respective basis states will be evenly populated.


Temperature in K, which is used with thermal population options SF_THERMAL and SO_THERMAL (keyword POPU). Default is set to 300 K.


Flag to include spin–orbit coupling. It is off by default.


Integer key governing the running mode of the program:

1 — start propagation from the RASSCF/RASSI HDF5 files (default);

2 — start propagation from reading intermediate file RDPREP. Can be used as a checkpointing/restarting option;

3 — only create the RDPREP without doing the propagation;


Time when propagation starts in fs. 0 fs by default


Time when propagation ends in fs. 10 fs by default.


Time step of integration in fs. Used only with fixed time step methods, see keyword METH. 0.0005 fs by default.


Method of integration: Runge–Kutta method of 4th order (classic_RK4) is set by default. RKCK (adaptive Runge–Kutta–Cash–Karp) with variable time step sometimes is better. Other available integrators are RK4, RK5, RK45 (adaptive Runge–Kutta–Fehlberg), and RK4_SPH (for propagation in spherical tensors basis).


Error threshold for the adaptive integration methods with variable step. \(10^{-6}\) by default.


Safety parameter to estimate error on each step for methods with adaptive step size. The estimated next step size is multiplied by the value of safety parameter. It is set to 0.9 by default. Increase it to get faster integration but at your own risk.


Time interval for writing populations and output of the discretized field. 0.05 fs by default.


Basis used for propagation, spin-free basis by default. For some features such as dipole moment or emission spectrum calculations SO basis is preferrable. Available options: CSF basis (could be convenient, e.g., for charge migration studies), SF or SO bases (could be useful to study spin dynamics). SPH enables the propagation in the basis of spherical tensors.


Basis used for the output of density matrix diagonal elements (populations), SF_SO means that density matrix is printed in two basis sets: spin-free and spin–orbit. Default value is set to SF. Available options: CSF, SF, SO, CSF_SF, SF_SO, CSF_SO, ALL.


Maximal rank of spherical tensor included in the spherical basis. Employed maximal rank can be reduced only if the SPH propagation basis is activated.


Maximal projection of spherical tensor included in the spherical basis. Employed maximal projection can be reduced only if the SPH propagation basis is activated.


Time step (in fs) for the output of the full density matrix. Without this keyword the full DM will NOT be saved in file RDOUT.


Flag to switch on the ionization. Works only if DYSOn keyword was used in RASSI.


Scaling parameter \(\alpha\) for the matrix of Dyson amplitudes.


A parameter defining the phenomenological decay rates of the ionized states.


Flag to switch on the dissipation due to the coupling to vibrational bath. See [137].


Number of vibrational modes included in the calculation. Needed only if IFDI is activated.


The coupling of the primary heat bath to the secondary one in \(\text{cm}^{-1}\). See [137], Eq. 7.


Enables reading of Huang–Rhys factors from file HR-FACT in spin–orbit coupled basis. See [137], Eq. 8. Provided that KEXT is not specified, both GAMM and HRSO are used to compute the dissipation rates according to Eq. 11 in [137].


Enables reading in the dissipation rate matrix \(k\) from the file KEXT. KEXTernal is an alternative to keywords HRSO and GAMM.


Number of incoming electric pulses, 1 by default. Set it to 0 if no pulse is needed.


Defines form of the pulse envelope function for each pulse. The electric field is supposed to be in the form \(A\vec{e}s(t)\sin{(\Omega(t-t_0)+\varphi_0)}\), where \(s(t)\) is the envelope function. Available options are

Gauss — Gaussian shape \(s(t)=\exp (-(t-t_0)^2/(2\sigma^2))\). Set by default.

sinX, cosX — more localized shape \(s(t)=\sin^X(\pi(t-t_0)/(2\sigma))\). Here X is the power to which the sine or cosine function is raised. For example, sin16 describes a \(\sin^{16}\) shape function. The support of a single pulse is \([t_0, t_0+2\sigma]\) for sinX and \([t_0-\sigma, t_0+\sigma]\) for cosX.

Mono — monochromatic pulse without shape function.

TYPE_X_CIRCLE — explicitly circularly polarized light, where X can be L or R for left and right direction, and TYPE replaces Mono or Gauss. For example, GAUSS_R_CIRCLE


On one line define as many amplitude values \(A\) in atomic units as many pulses you ask for.


Here should be shifts \(t_0\) in fs for each pulse with respect to the global initial time point.


Three complex numbers in the format (XR,XI) (YR,YI) (ZR,ZI) defining the polarization \(\vec{e}\). By default, the electric field is considered to be linear polarized along the \(x\)-direction. If the number of pulses is more than one, the polarization vector should be given for each pulse on a separate line.


Pulse dispersion \(\sigma\) in fs for each pulse (in one line separated by space). See keyword PTYP for definition.


Carrier frequency \(\Omega\) in eV for each pulse (in one line separated by space).


Phase \(\varphi_0\) in radians for each pulse.


Enable correction to carrier frequency simulating experimental linear frequency chirp. Due to time-dependent phase, the carrier frequency gets an additional linear term \(\Omega \rightarrow \Omega + a (t - t_0)\). The constant \(a\) should be specified.


Enable correction to electromagnetic pulse as if is given as time derivative of a vector potential.


Activates the calculation of mean value of dipole moment, currently it is written to the file DIPOLE.dat


Activates the calculation of emission spectrum Input examples

> copy /path/to/file/s3.rasscf.h5  RASSD1
> copy /path/to/file/s1.rasscf.h5  RASSD2
> copy /path/to/file/si.rassisd.h5 RASSISD
> copy /path/to/file/kmatrix.dat   HRFACT


NRSManifolds          = 2
NRDEt,CSF,STATES,SPIN = 25   25   25   3
                        30   30   30   1
NSTAte                = 105 all
FINAltime             = 10
AMPLitude             = 9.0
TAUShift              = 1.
SIGMa                 = 0.3
OMEGa                 = 875

NRSManifolds          = 2
NRDEt,CSF,STATES,SPIN = 25   25   25   3
                        30   30   30   1
POPUlatedstyle        = SO_THERMAL
NSTAte                = 105 all
FINAltime             = 6
Tout                  = 0.0005
METHod                = RKCK
DMBAsis               = SO
PROPbasis             = SO
TFDM                  = 0.005

PTYPe                 = Gaussian
NPULses               = 1
AMPLitude             = 9.0
TAUShift              = 1.
POLArization          = (1.0,0.0) (0.0,0.0) (0.0,0.0)
SIGMa                 = 0.3
OMEGa                 = 875
PHASe                 = 0