TURBOMOLE manual
v8.0

7.3 How to Perform a Calculation

7.3.1 Basis Sets for Periodic Calculations

In periodic systems, overly diffuse Gaussian basis functions can lead to near-singular overlap matrices and unstable SCF convergence. As a rule of thumb, primitives with exponents smaller than 0.01 should be avoided, and exponents larger than 0.1 are preferable. A straightforward remedy is to remove the most diffuse primitives from standard molecular basis sets.

A more robust approach is to employ basis sets specifically optimized for periodic DFT calculations. TURBOMOLE provides the pob-DZVP-rev2 and pob-TZVP-rev2 families, which are designed for solid-state calculations and offer improved stability and accuracy compared to unmodified molecular basis sets [209, 210, 211, 212]. These basis sets are available in the TURBOMOLE basis library and are recommended for all periodic calculations.

7.3.2 Prerequisites

riper calculations require a properly formatted control file containing essential system information and calculation parameters.

Initial Density and Charge Specification

By default, riper constructs the initial electron density by superposing atomic densities (see Sec. 25.2.2) and assumes a charge-neutral system. This superposition provides a robust and transferable starting point and is recommended for most applications. The parameters can be customized in the $atomdens group to specify:

  • Total system charge

  • Spin multiplicity for open-shell systems

  • Custom atomic charges if needed

Alternatively, define can be used to generate initial orbitals and define both the total charge and orbital occupations explicitly (see Sec. 4).

Basis Set Specification

Both the orbital basis set and the auxiliary basis set for density fitting must be specified in the $atoms group. The universal auxiliary basis set is typically sufficient for most applications.

Minimal Control File Examples

A minimal control file for a neutral water molecule without explicit charge specification:

	$atoms
	  basis =pob-DZVP-rev2
	  jbas  =universal
	$coord
	  -0.12095422215797   0.00000000000000   0.00000000000000   o
	   1.91995057443318   0.00000000000000   0.00000000000000   h
	  -0.80125579971337   0.00000000000000   1.92418350276796   h
	$end

For a charged system (+1) with open-shell configuration:

	$atoms
	  basis =pob-DZVP-rev2
	  jbas  =universal
	$coord
	  -0.12095422215797   0.00000000000000   0.00000000000000   o
	   1.91995057443318   0.00000000000000   0.00000000000000   h
	  -0.80125579971337   0.00000000000000   1.92418350276796   h
	$atomdens
	  charge 1
	$uhf
	$end

Note: If $uhf is not specified for an odd number of electrons, riper assumes a closed-shell configuration with fractional occupation.

Exchange-Correlation Functional Selection

By default, a Hartree–Fock single-point energy calculation is performed. For DFT calculations, specify the desired functional in the $dft data group. riper supports:

  • Local density approximation (LDA) functionals

  • Generalized gradient approximation (GGA) functionals

  • Meta-GGA functionals

  • Global and range-separated hybrid functionals

  • Functionals available via the LibXC and XCFun libraries (see Sec. 6.3)

For periodic systems, additional keywords specifying periodicity and \(k\)-point sampling must be added (detailed in the following sections).

7.3.3 Creating the Input File

This section provides comprehensive guidance for creating riper input files for systems ranging from molecules (0D) to fully periodic crystals (3D). While the complete control file must be created manually or using define, the RIPER-Tools web interface can assist with structure preparation, k-point generation, and output analysis.

Using RIPER-Tools for Structure Preparation and Input Assistance

Navigate to https://ripertools.turbomole.org/ and use it for:

  1. Import and convert structures: Upload structure files (CIF, POSCAR, XYZ, Quantum ESPRESSO) or retrieve from Materials Project/PubChem databases, with automatic conversion to TURBOMOLE’s coord format

  2. Manipulate structures: Create supercells, translate structures within unit cells, place adsorbates on surfaces, or pack molecules in periodic cells

  3. Generate \(k\)-point paths: Automatically determine high-symmetry \(k\)-points for band structure calculations based on the Bravais lattice, generating the $kpoints snippet

  4. Create RT-TDDFT inputs: Generate RT-TDDFT input files with appropriate electric field parameters

  5. Analyze output: Parse and visualize RT-TDDFT spectra, band structures, density of states, and other output files

RIPER-Tools significantly reduces preparation time by handling format conversions, automatically generating k-point paths for band structures, and providing various structure manipulation capabilities. The generated coord files and $kpoints snippets can be directly incorporated into manually created or define-generated control files.

Manual Input File Creation

For users requiring fine control over calculation parameters, manual input file creation follows these steps:

  1. Generate the control file: Create the control file either:

    • Manually, following the prerequisites in Sec. 7.3.2

    • Using define (see Sec. 4) — ensure \(C_1\) symmetry is selected in the geometry menu, as riper currently supports only \(C_1\) symmetry

  2. Specify system periodicity: riper defaults to molecular (0D) calculations. For periodic systems, add:

    • $periodic \(n\) where \(n = 1, 2, 3\) for 1D, 2D, or 3D periodicity

    • Unit cell parameters via $cell or lattice vectors via $lattice

  3. Define k-point sampling: For periodic calculations, specify the k-point grid using $kpoints. RIPER-Tools can assist by generating k-point paths for band structure calculations.

  4. Configure optional parameters: Add the $riper data group with calculation-specific keywords (see Sec. 25.2.15)

Note: Relativistic calculations with effective core potentials require no additional keywords compared to ridft (see Sec. 6.5).

Unit Cell Specification

Unit cells can be defined using either conventional crystallographic parameters or explicit lattice vectors. By default, lengths are in atomic units (bohr) and angles in degrees. To use angstroms, append angs to the keyword (e.g., $cell angs or $lattice angs).

Method 1: Cell Parameters ($cell)

The required parameters depend on system dimensionality:

3D systems require six parameters:

  • \(|\mathbf{a}|\), \(|\mathbf{b}|\), \(|\mathbf{c}|\): lengths of cell vectors

  • \(\alpha\): angle between vectors \(\mathbf{b}\) and \(\mathbf{c}\)

  • \(\beta\): angle between vectors \(\mathbf{a}\) and \(\mathbf{c}\)

  • \(\gamma\): angle between vectors \(\mathbf{a}\) and \(\mathbf{b}\)

  • Convention: \(\mathbf{a}\) aligns with the \(x\)-axis, \(\mathbf{b}\) lies in the \(xy\)-plane

2D systems require three parameters:

  • \(|\mathbf{a}|\), \(|\mathbf{b}|\): lengths of cell vectors

  • \(\gamma\): angle between \(\mathbf{a}\) and \(\mathbf{b}\)

  • Convention: \(\mathbf{a}\) aligns with the \(x\)-axis, \(\mathbf{b}\) lies in the \(xy\)-plane

1D systems require one parameter:

  • \(|\mathbf{a}|\): length along the \(x\)-axis

Method 2: Lattice Vectors ($lattice)

Alternatively, provide explicit lattice vectors:

  • 3D systems: Three 3D vectors (nine components)

  • 2D systems: Two 2D vectors in the \(xy\)-plane (four components)

  • 1D systems: One scalar value (length along \(x\)-axis)

TURBOMOLE also comes with the cif2t script that can be used to create the coord file for periodic systems from a CIF file. This can be used as:

	cif2t my_cell.cif

Basic Examples

Example 1: 1D Periodic System

Unit cell with \(|\mathbf{a}| = 5\) bohr:

	$periodic 1
	$cell
	  5.0

Equivalent using lattice vectors:

	$periodic 1
	$lattice
	  5.000000000000000

Example 2: 2D Periodic System

Surface cell with \(|\mathbf{a}| = 5\) bohr, \(|\mathbf{b}| = 8\) bohr, \(\gamma = 60^o\):

	$periodic 2
	$cell
	  5.0 8.0 60.0

Equivalent using lattice vectors:

	$periodic 2
	$lattice
	  5.000000000000000  0.000000000000000
	  6.928203230275509  4.000000000000000

Example 3: 3D Triclinic System

Unit cell with \(|\mathbf{a}| = 5\), \(|\mathbf{b}| = 8\), \(|\mathbf{c}| = 6\) bohr; \(\alpha = 70^o\), \(\beta = 60^o\), \(\gamma = 90^o\):

	$periodic 3
	$cell
	  5.0 8.0 6.0 70.0 60.0 90.0

Equivalent using lattice vectors:

	$periodic 3
	$lattice
	  5.000000000000000  0.000000000000000  0.000000000000000
	  0.000000000000000  8.000000000000000  0.000000000000000
	  3.000000000000000  2.052380899954828  5.637041127056543

Complete Input Examples

Example 1: 3D Cubic Crystal with K-Point Sampling

This example demonstrates a periodic bulk calculation with metallic character:

	$coord file=coord
	$periodic 3
	$atoms
	  basis =pob-TZVP-rev2
	  jbas  =universal
	$cell
	  18.5911 16.5747 16.5747 90. 90. 90.
	$dft
	  functional pbe
	$kpoints
	  nkpoints 3 3 3
	$riper
	  lenonly on    # Energy-only calculation (skip gradients)
	  sigma 0.01    # Gaussian smearing for metallic systems
	$end

Key features:

  • Cubic unit cell with PBC in all directions

  • \(3 \times 3 \times 3\) k-point grid for BZ sampling

  • PBE functional for solid-state DFT

  • Gaussian smearing to handle metallic occupation

Example 2: 2D Surface with Enhanced Multipole Accuracy

Surface calculation with modified CFMM parameters:

	$coord file=coord
	$periodic 2
	$atoms
	  basis =pob-TZVP-rev2
	  jbas  =universal
	$cell angs
	  18.5911 16.5747 90.0
	$dft
	  functional pbe
	$kpoints
	  nkpoints 3 3
	$riper
	  lmaxmom 30    # Increased multipole order (default: 20)
	$end

Key features:

  • 2D periodic slab geometry

  • Cell parameters in angstroms

  • Enhanced multipole expansion for improved accuracy

Example 3: 1D Nanowire

One-dimensional periodic system:

	$coord file=coord
	$periodic 1
	$atoms
	  basis =pob-DZVP-rev2
	  jbas  =universal
	$cell
	  18.5911
	$dft
	  functional pbe
	$kpoints
	  nkpoints 3
	$end

Example 4: Large Molecular Complex

Non-periodic calculation with memory optimization:

	$coord file=coord
	$atoms
	  basis =pob-DZVP-rev2
	  jbas  =universal
	$dft
	  functional pbe
	$riper
	  lpcg on    # Low-memory PCG algorithm
	$end

Example 5: Custom \(K\)-Point Sampling

User-defined \(k\)-points for band structure analysis:

	$coord file=coord
	$periodic 3
	$atoms
	  basis =pob-TZVP-rev2
	  jbas  =universal
	$cell
	  16.5747 16.5747 16.5747 90. 90. 90.
	$dft
	  functional pbe
	$kpoints
	  custom 5
	  # k_x    k_y    k_z   weight
	  0.00   0.00   0.00   1    # Gamma point
	  0.25   0.25   0.25   8    # L point
	  0.50   0.50   0.50   1    # R point
	  0.00   0.50   0.00   6    # X point
	  0.50   0.00   0.50   3    # M point
	$end

Important Notes:

  • \(k\)-point weights are internally normalized to sum to 1.0

  • Consult crystallographic references for optimal k-point sets for specific Bravais lattices

  • For metals, always use smearing techniques (see Sec. 7.4.4)

  • Universal auxiliary basis is typically sufficient but can be customized if needed

  • All examples assume structures are provided in separate coord files

  • RIPER-Tools can automatically generate these input files through its web interface

7.3.4 Single Point Energy and Gradient

By default riper calculates single point energy and gradient in one run. For this, simply invoke riper as

	nohup riper > riper.out &

For energy calculation only add lenonly on to the $riper section in the control file, i.e.,

	$riper
	  lenonly on

The parallel OpenMP version of riper can be invoked as described in Sec. 3.4.2.

7.3.5 Structure Optimization

jobex will automatically use riper when the keyword $periodic is present in the control file. Alternatively, the use of riper can be forced by specifying -riper argument of jobex, i.e., invoking

	nohup jobex -riper > jobex.out &

Note, that at present only structure optimization using Cartesian coordinates is possible when using periodic boundary conditions. Using internal coordinates for periodic systems will lead to convergence problems in structure optimizations.

7.3.6 Optimization of Cell Parameters

Simultaneous optimization of atomic positions and cell parameters or lattice vectors can be performed by specifying the keyword $optcell in the control file and then invoking the jobex script as described in 7.3.5.

The calculation of energy first derivatives with respect to lattice vectors usually requires higher accuracy than nuclear gradients. Therefore, in case of convergence issues it is recommended to decrease the SCF convergence threshold to at least \(1.0\times 10^{-7}\) by specifying $scfconv 7. In addition, convergence can often be improved by including weight derivatives by adding the option weight derivatives in the $dft data group.

7.3.7 Band Structure Plots

Band structure plots show the values of band (orbital) energies for values of \(\textbf{k}\) along lines connecting symmetry points. In riper these lines can be specified by providing their start and end points as well as the number of \(\textbf{k}\) points along the line within the $kpoints section. They will be calculated at the end of SCF procedure and written to the file bands.xyz.

In the following example band energies are calculated along four lines, as specified by the keyword kptlines 4. Each line definition starts in a new line with the keyword recipr, followed by three real numbers defining the start point of the line and three real numbers defining its end point. Finally, the number of \(\textbf{k}\) points along the line is given as an integer number. Thus, the first line starts at the point (0.500 0.500 0.500), ends at (0.000 0.000 0.000) and contains 40 \(\textbf{k}\) points.

	$kpoints
	  kptlines 4
	  recipr 0.500 0.500 0.500 0.000 0.000 0.000 40
	  recipr 0.000 0.000 0.000 0.500 0.500 0.000 40
	  recipr 0.500 0.500 0.000 0.746 0.373 0.373 40
	  recipr 0.746 0.373 0.373 0.000 0.000 0.000 40

Each line of the resulting bands.xyz file contains five real numbers: the coordinates \(k_1\), \(k_2\) and \(k_3\) of the \(\textbf{k}\) point, its length \(|\textbf{k}|\) and the corresponding band energy \(\epsilon_{n\mathbf{k}}\).

TURBOMOLE comes with the kpath_coord script that can help with generating the kpath for band structure calculations, in addition to RIPER Tools. Given a coord file with atomic positions in atomic units (Bohrs) and cell parameters either in the coord file or the control file, the script can be invoked as follows

	kpath_coord coord

Running the above command will give the snippet required to be added to the control file for band structure calculations. Note, for the band structure calculation of pristine materials, one should use the primitive unit cell to avoid band folding effects. The above script can also help with the conversion of a conventional unit cell to the primitive unit cell.

7.3.8 Calculation of Densities and MOs on Grids

How to Perform

Calculation of data on grids for molecular and periodic systems using riper requires the keyword $pointvalper in the control file. The values are calculated on a 3D grid and written to appropriate output files. These files can be generated either running a single point energy riper calculation or invoking

	nohup riper -proper > riper.out &

provided that converged orbitals and occupation numbers from previous calculation are present in the files RIPER.BANDS and RIPER.BANDS.OCCUPATIONS, respectively.

The format of the output files can be specified using the keyword fmt= following $pointvalper in the same line. Supported formats are:

plt

(default) The generated data is written to binary files that can be read by gOpenMol and other external visualization programs. This format uses orthogonal grids. Therefore, for non-orthogonal unit cells grid data is generated for a rectangular box that contains the supercell (unit cell and its periodic images). By default, the values at grid points outside of the supercell are set to zero. For strongly non-orthogonal systems this may lead to large files. The option full switches off the zeroing of values on grid points outside the supercell.

cub

Gaussian cube format. Can be imported by several external visualization programs.

xyz

Coordinates of grid points and calculated values are stored in a text file. Each line contains Cartesian grid point coordinates and the corresponding value.

upt

Values and grid data is output to binary files. First 8 (integer) records are: 1-2) rank value and descriptor, 3-8) number of grid points and number of periodic unit cell images in each periodic direction, 9...) Cartesian coordinates and the corresponding value repeated for all grid points.

Output files in cub format contain information about atomic coordinates. For other formats additional coord.xyz file containing atomic coordinates is generated.

The following options can be used along with the keyword $pointvalper:

nimg n1 n2 n3

Number of unit cell images n1, n2 and n3 in the periodic directions a, b and c, respectively, for which plot data is generated.

npts n1 n2 n3

Number of grid points n1, n2 and n3 along each periodic direction. If not specified, value 100 is used for each n.

eps real

Specifies the distance real in bohr around the system for which plot grid is generated in aperiodic directions. Default value is 5 bohr.

ngrdpbx n

Number of grid points stored in one octree box during density calculations. Default value is 50. For very large systems or high resolution it may be necessary to increase this parameter to avoid memory allocation problems.

full

Only valid for plt output format. Switches off zeroing of values on grid points outside the supercell.

Electron Density

Plots of electron density require option dens in the $pointvalper data group

$pointvalper 
  dens

The values of the total density and, in case of UHF calculations also the spin density are written to files td.fmt and sd.fmt, respectively, where fmt is the selected format.

Molecular Orbitals

Plot files for visualization of molecular orbitals can be generated by setting orbs n in the $pointvalper data group, where n is the number of orbitals to plot. The way to specify the orbitals is best explained by some examples:

example 1:

An open shell system is considered. Number of \(\mathbf{k}\) points along the reciprocal lattice vectors \(\mathbf{b}_1\), \(\mathbf{b}_2\) and \(\mathbf{b}_3\) is 3, 5, and 7, respectively. Possible \(\mathbf{k}\) point components \(k_j\) calculated from equation 7.20 are therefore \(k_1 \in \left( \frac{-1}{3}, 0, \frac{1}{3} \right)\), \(k_2 \in \left( \frac{-2}{5}, \frac{-1}{5}, 0, \frac{1}{5}, \frac{2}{5} \right)\), and \(k_2 \in \left( \frac{-3}{7}, ... ,\frac{3}{7} \right)\) The kpoints and $pointvalper groups are:

$kpoints
  nkpoints 3 5 7
$pointvalper 
  orbs 2
  k 1 4 6 a 1 
  k 2 3 4 b 3

Here, two orbitals will be plotted. They are defined in two lines following the orb 2 option. Each line starts with k followed by three numbers that specify which of possible components \(k_j\) are used. In this example \(\mathbf{k}\) point components for the first orbital are \(\left( \frac{-1}{3}, \frac{1}{5}, \frac{2}{7} \right)\). For the second one \(\mathbf{k}=(0,0,0)\), i.e., this orbital is at the \(\Gamma\)-point. Subsequently, type and number of the orbitals are defined (orbitals at each \(\mathbf{k}\)-point are ordered by increasing energy). Here the first one is a 1 - the lowest \(\alpha\) orbital. Similarly, b 3 denotes the third lowest energy \(\beta\) orbital. For a non-\(\Gamma\) \(\mathbf{k}\) point, the calculation will generate three plottable files corresponding to the: (i) real component of \(\Psi\), (ii) imaginary component of \(\Psi\), and (iii)\(\Psi^2\), for a given orbital. Note: Only the real components are plotted for the \(\Gamma\)-point.

example 2: In this example a closed shell system with 3, 2 and 4 \(\mathbf{k}\) points in each direction is given:

$kpoints
  nkpoints 3 2 4
$pointvalper 
  orbs 3
  k 1 2 1 a 1
  k 0 0 0 a 2

Possible \(\mathbf{k}\) point components from equation 7.20 are \(k_1 \in \left(\frac{-1}{3},0,\frac{1}{3}\right)\), \(k_2 \in \left( \frac{-1}{4},\frac{1}{4} \right)\) and \(k_2 \in \left( \frac{-3}{8},\frac{-1}{8},\frac{1}{8},\frac{3}{8} \right)\). In analogy to the previous example, three data files for the real component, imaginary component, and square of the wavefunction (MO) will be generated for the lowest energy orbital at \(\mathbf{k}\) point \(\left( \frac{-1}{3},\frac{1}{4},\frac{-3}{8} \right)\). For the second molecular orbital only the real component will be generated at the \(\Gamma\)-point. Orbitals at the \(\Gamma\)-point are defined by using ‘special’ 0 0 0 component indices, as it is not included in the grid for even number of \(\mathbf{k}\) points.

example 3: In this example an input for a \(\Gamma\)-point calculation is shown:

$pointvalper fmt=plt
  orbs 2  
  k 1 1 1 a 5 
  k 1 1 1 b 5

Here k 1 1 1 is the only valid option. For \(\Gamma\)-point calculation, orbitals are real, thus it is the only component that is plotted. In this case the fifth lowest energy orbitals with spin components \(\alpha\) and \(\beta\) are plotted. The output format is set to plt using the fmt= keyword.

example 4: In this example an input for a custom \(\mathbf{k}\) point calculation is shown:

$kpoints
  custom 2
  0.00  0.00  0.00   1
  0.50  0.50  0.25   1
$pointvalper fmt=cub
  orbs 2
  k 1 a 5 
  k 2 a 6

Instead of three indices after k (like in the previous examples), only a single index is required, specifying the index of the custom \(\mathbf{k}\) point in the list given by the user. The above example will plot the fifth lowest energy molecular orbital at the first \(\mathbf{k}\) point (\(\Gamma\)) and the sixth MO at the second, (0.50 0.50 0.25) \(\mathbf{k}\) point. Note: If the list of custom \(\mathbf{k}\) points doesn’t include a \(\Gamma\)-point, one can still plot the MOs at this point, as a \(\Gamma\)-point with 0 weight is prepended to the supplied list with an index 0. In such a case, to plot the fifth lowest energy orbital with spin \(\alpha\) at the \(\Gamma\)-point, the input would be k 0 a 5 .

7.3.9 Density of States

A simulated density of states (DOS) is calculated by broadening the discrete energy levels with Gaussians and superimposing them. The output files (dos for RKS calculations, and dos_a+b, dos_a-b, dos_alpha, dos_beta in case of UKS) contain energies (first column) and corresponding total DOS (second column) as well as s-, p-, ... contributions (following columns). The calculation is controlled by the keyword $dosper:

$dosper width=real emin=real emax=real scal=real npt=integer

The parameters in example above have the following meaning:

the width of each Gaussian, default value is 0.01 a.u.

lower/upper bounds for energy in DOS calculation

scaling factor for DOS (total and s-, p-, ... contributions)

resolution (number of points)

All parameters are optional, default values are usually sufficient to generate a satisfactory plot. DOS files can be generated either running a single point energy riper calculation or using -proper option

nohup riper -proper > riper.out &

provided that converged orbitals and occupation numbers from previous calculation are present in the files RIPER.BANDS and RIPER.BANDS.OCCUPATIONS, respectively.

7.3.10 RT-TDDFT

How to Perform

RT-TDDFT and RT-TDHF calculations require data-groups: $fields, $electric field and $rttddft to provide the information about the type of field, electric field and the RT-TDDFT calculation, respectively. Additionally, keywords $rtdipol and $rtenergy followed by file names may be provided in the control file to print out the dipole moment and the energies respectively every \(n\) number of steps. It is recommended to monitor the values in these files as the calculation is running to ensure that the simulation is stable and not diverging. To calculate the absorption spectrum and print it, the keyword $rtspectrum followed by the units is needed. Possible values of the units are ev, nm, cm-1, and hartree. Example:

$rtspectrum ev

For now, only electric fields are supported. So, the $fields block would look like

$fields
  electric on

The above would enable the use of an electric field.

The information pertaining to the electric field is provided within the $electric field data-group.

There are currently three types of electric fields supported.

  • Gaussian

  • Static

  • Laser

Static fields only require the Cartesian components of the amplitude. The Gaussian field requires two more parameters: peak position and peak width. Finally, the Laser field requires the Cartesian components, frequency and the full width at half maximum (FWHM) of the field envelope. Additionally, the phase may also be provided.

Note: The absorption spectrum can only be calculated for the Gaussian field.

To define a static field

$electric field
  amplitude x=2.0E-5 y=2.0E-5 z=2.0E-5
  static

To define a Gaussian ("kick") field

$electric field
  amplitude x=2.0E-5 y=2.0E-5 z=2.0E-5
  gaussian tzero=3.0 width=0.2

To define a Laser field

$electric field
  amplitude x=0.005d0 y=0.0050 z=0.0d0
  phase x=1.570796d0 y=0.0d0 z=0.0d0
  laser omega=0.056961 sigma=750.0d0

Note: the amplitudes are provided in atomic units (1 au \(=\) 514.2 V/nm). The Gaussian electric field is defined as \[\begin{equation} \vec{E}(t)=\vec{A} \mathrm{e}^{-\frac{\left(t-t_{0}\right)^{2}}{2 w^{2}}} \end{equation}\](7.1)

tzero in the above code snippet corresponds to \(t_0\) (center of the peak) and width corresponds to \(w\) (which gives a measure of the width of the pulse). Note: both are in atomic units of time (1 au \(=\) 0.02419 fs)

The Laser electric field is defined as \[\begin{equation} \mathbf{E}(t)=f(t) \left(E_{x} \sin \left(\omega_{0} t+\phi_x\right)\mathbf{n_x}+E_{y} \sin \left(\omega_{0} t+\phi_y\right)\mathbf{n_y}+E_{z} \sin \left(\omega_{0} t+\phi_z\right)\mathbf{n_z}\right) \end{equation}\](7.2) where \(E_{x,y,z}\) are the amplitudes along \(x,y,z\) axes, \(\omega_0\) is the carrier frequency, \(\phi_{x,y,z}\) are the carrier-envelope phases, and \(f(t)\) is an envelope/shape function given as \[\begin{equation} f(t)= \begin{cases}\cos ^{2}\left(\frac{\pi}{2 \sigma}(t-\sigma)\right) & \text { if } 0 \leq t \leq 2 \sigma \\ 0 & \text { otherwise }\end{cases} \enspace , \end{equation}\](7.3) where \(\sigma\) is the FWHM of the field envelope. The total pulse duration is \(2\sigma\). The laser pulse reaches its maximum at \(t=\sigma\), corresponding to a mean intensity of \[\begin{equation} I=\frac{1}{2} \varepsilon_{0} c | E |^{2} \enspace (\text{SI units}), \end{equation}\](7.4) where \(\varepsilon\) is the permittivity of vacuum, \(| E |^{2}=E_x^2+E_y^2+E_z^2\), and \(c\) is the velocity of light.

omega in the above code snippet for Laser field corresponds to \(\omega\) the frequency of the laser pulse and sigma corresponds to the FWHM of the laser envelope. The expected units for frequency is a.u. (Note: 1 eV \(= 0.0367493\) a.u.). Similarly, sigma is expected in a.u. of time (1 au \(=\) 0.02419 fs). Finally, phase corresponds to the phase angle \(\phi\) in eq. 7.2. It is expected in radians (1 rad\(= 57.2958^o\)). Phase is useful to create a circularly polarized laser pulse.

Finally, the parameters regarding the RT-TDDFT calculation go within the $rttddft data-group. The meanings are explained below

magnus

Can take values 2 or 4. "2" for second order Magnus expansion and "4" for fourth order Magnus expansion. Default value is 2.

scf

if on, then SCF procedure is used for the time integration. If off
then Predictor-Corrector scheme is used instead.

iterlim

Max SCF cycles if scf is on. Default value is 15.

time

Specifies the evolution time in au. (1 au\(=\) 0.02419 fs)

tstep

The time step for the time evolution in au. 0.1 au is usually a good starting point.

print step

Specifies the number of steps \(n\) after which the dipole moments and energies are printed out if requested. Default value is 100. That means the quantities are printed out at every 100 steps. To have all the information for post processing, a value of 1 is recommended.

damping

Only valid for absorption spectrum calculation. It is the factor gamma in the equation to calculate the complex polarizability tensor. Default value is 0.004 au. Recommended values in the range of 0.003 au to 0.005 au.

min energy

Only valid for absorption spectrum calculation. Specifies the minimum value of the energy range used to perform the Fourier transform from time to frequency space. Units: au. Default value is 0.15 au.

max energy

Only valid for absorption spectrum calculation. Specifies the maximum value of the energy range used to perform the Fourier transform from time to frequency space. Units: au. Default value is 0.625 au.

energy step

Only valid for absorption spectrum calculation. Specifies the step value or energy interval dE at which to sample the energy values for Fourier transform and absorption spectrum plotting. Units: au. Default value is 0.005 au.

Consider the following example.

Example 1:

Here, an absorption spectrum calculation for a molecule is performed using RT-TDDFT. The following keywords are expected in the control file.

$fields
  electric on
$electric field
  amplitude x=2.0E-5 y=2.0E-5 z=2.0E-5
  gaussian tzero=3.0 width=0.2
$rttddft 
  magnus 2
  scf off
  time 1000.0d0
  tstep 0.1d0
  print step 1 
  min energy = 0.013d0 
  max energy = 0.5d0
  energy step 0.001d0
$rtdipol
$rtenergy
$rtspectrum ev

The data-group $fields is used to specify that only electric field is needed using electric on. Here, a Gaussian "kick" pulse is used to excite the system. The parameters for which are specified in the $electric field section. The pulse has a small amplitude in all three directions. The amplitude is specified in au (1 au \(=\) 514.2 V/m). The position of the maximum of the pulse is at \(t=\)3 au, specified using tzero. The width of the pulse is specified using width parameter Eq. (7.1). The simulation goes on for 1000 au (1 au\(=\) 0.02419 fs) in steps of 0.1 au as specified by the time and tstep keywords in the $rttddft section, respectively. Second order Magnus expansion is used for the time-evolution operator indicated by magnus 2. Predictor-corrector scheme is used for time-integration since scf is off. Notice here, that we are also requesting to print out energies and dipole moments using the $rtenergy and $rtdipol keywords, respectively. The default names of the corresponding files are rtenrgy and rtdipo. To print out to a custom file name use $rtenergy file filename instead. It is recommended to have these files printed as they allow to monitor the energies in real time and an unstable simulation can be detected early. The print step 1 in the $rttddft section specifies that the dipole moments and the energies are printed at each iteration, i.e. every step. This is useful if post-processing is needed after the calculation. Default value of print step is 100. Finally, to calculate the spectrum, we provide the $rtspectrum keyword followed by the units (in this case electron volts). The energy range for which the absorption spectrum is calculated is 0.35 eV (0.013 au) to 13.6 eV (0.5 au), as specified using the min energy and max energy keywords with the step interval of 0.001 au. The spectrum is printed out in a file called rtspec.
The calculation is performed by running
nohup riper > riper.out &
To keep a track of the progress of the calculation, i.e. how many steps have been completed, check the rtenrgy file.
The RT-TDDFT absorption spectrum can be plotted from the rtspec file after the calculation ends.

Example 2:

In this example the electron density is plotted at each time step, to visualize the electron dynamics in real-time. The following keywords are expected in the control file.

$fields
  electric on
$electric field
  amplitude x=2.0E-5 y=2.0E-5 z=2.0E-5
  gaussian tzero=3.0 width=0.2
$rttddft
  magnus 2
  scf off
  time 1000.0d0
  tstep 0.5d0
  print step 1
$rtdipol
$rtenergy
$rtdens
$pointvalper fmt=cub
  dens

Two new keywords are noticed here compared to example 1. $rtdens keyword tells the program to print out the difference of the excited state and ground state density. Second, the $pointvalper section, explained already, is used to specify the format in which the electron density is printed out.

7.3.11 DFT-based Embedding

How to Perform

In order ro perform DFT-based embedding using the riperembed python script one requires two files: (i) an input file with the name input and (ii) a coordinate file with the name totalCoord with the atomic coordinates of the total system in TURBOMOLE format (bohrs). The order of atomic coordinates should be such that the atoms of the active subsystem are in the beginning followed by the atoms of the environment subsystem.

To facilitate the creation of the input file a web application is made available at this url: https://ripertools.turbomole.org/ The supported options for the input file are provided below

nsystm

(integer)
Number of systems (Currently only 2 subsystems are supported)

ptnIndx

(2 element list of integers)
Atomic indices of the last atom of the two subsystems (Ex: for a system of HF dimer with one HF being the active subsystem and the other monomer being the environment, the input value would be 2 4)

KEfunc

(integer or string)
LibXC code of the KEDF (Ex: 521 for LC94/PW91k) or if one doesn’t want to use a KEDF then one can set the value to electro as then the embedding can be thought of electrostatic embedding

basis

(string)
Basis set for the total system or environment. (Default: def2-SVP)

basis_cluster

(string; optional)
Basis set for the active subsystem. If this option is not specified then the basis set specified in basis is used for both active and environment subsystems (Ex: def2-TZVP)

auxbasis

(string; optional)
Auxiliary basis set for both active and environment subsystems (Default: universal)

xName

(integer)
LibXC code of the Exchange functional (Default: 101 for GGA_PBE_X)

cName

(integer)
LibXC code of the Correlation functional (Default: 130 for GGA_PBE_C)

frozen

(integer)
Perform Freeze-and-Thaw (FaT) if True (Default: False)

super

(string)
Use supermolecular basis if True (Default: False)

emb_error

(string)
Calculate error in total energy with respect to a KS-DFT calculation on the total system. Compatible with method=1,3,5. (Default: True)

method

(integer)
Specify the embedding method. Supported options are 1 for molecule-in-molecule FDE, 2 for molecule-in-periodic FDE, 3 for molecule-in-molecule PbE, 4 for embedding using a fixed embedding potential read from the file embdpot, 5 for periodic-in-periodic PbE (only supermolecular basis is supported for method=5)

scfconv

(integer)
SCF convergence threshold (Default: 7)

denconv

(integer)
Density convergence threshold (Default: 7)

ricore

(integer)
Memory size for ricore (Default: 2000)

nmax_FaT

(integer)
Maximum number of FaT cycles for FDE and PbE (Default: 5)

mxitdiis

(integer)
DIIS subspace (Default: 5)

scfiterlimit

(integer)
SCF max iteration limit (Default: 500)

charge_cluster

(integer)
Active subsystem charge, but make sure that the subsystem is closed-shell (Default: 0)

charge_env

(integer)
Environment subsystem charge, but make sure that the subsystem is closed-shell (Default: 0)

cell_params

(list of float)
For periodic-in-periodic or molecule-in-perriodic calculations, the cell parameters need be specified in Angstroms and degrees (Ex: 5.0 5.0 5.0 90.0 90.0 90.0 for 3D PBC, 5.0 5.0 90.0 for 2D PBC, 5.0 for 1D PBC)

periodicity

(integer)
Periodicity for periodic calculations, which can be 1, 2, 3 (Ex: 3)

kpoints

(string + list of integers)
Number of \(k\)-points for periodic calculations (Ex: nkpoints 3 3 3)

The embedding calculation is run as follows

nohup riperembed > outputEmbed &

The following examples illustrate specific-use cases:

Example 1 FDE for molecule-in-molecule embedding (monomolecular basis):

Consider the system of an HF dimer. One HF molecule will serve as active subsystem and the other as environment. The totalCoord file looks like

$coord
   2.23021697819493     -0.20903016559004     -0.00001133835675      f
   3.01766585447976      1.34738606547718      0.00002645616575      h
  -2.75370890934066      0.03639612516738     -0.00001889726125      f
  -1.01914819647029     -0.17888147499189      0.00028345891875      h
$end

The coordinates of the active subsystem are provided in the beginning. In order to perform an FDE calculation for molecular active and environment subsystems, the input file would look like the following

$EMBED
  nsystm = 2
  ptnIndx = 2 4
  #LC94
  KEfunc = 521
  basis = def2-SVP
  auxbasis = universal
  #PBE XC
  xName = 101
  cName = 130
  frozen = True
  super = False
  emb_error = True 
  method = 1
  ricore = 2000
  scfiterlimit = 30

Example 2 PbE for molecule-in-molecule embedding (supermolecular basis):

In order to perform an PbE calculation for molecular active and environment subsystems using a supermolecular basis, the input file would look like the following

$EMBED
  nsystm = 2
  ptnIndx = 2 4
  #No KEDF
  KEfunc = electro
  basis = def2-SVP
  auxbasis = universal
  #PBE XC
  xName = 101
  cName = 130
  frozen = True
  super = True
  emb_error = True 
  method = 3
  ricore = 2000
  scfiterlimit = 30

Example 3 FDE+FaT for molecule-in-molecule embedding (monomolecular basis):

In order to perform an FDE+FaT calculation for molecular active and environment subsystems, the input file would look like the following

$EMBED
  nsystm = 2
  ptnIndx = 2 4
  #LC94
  KEfunc = 521
  basis = def2-SVP
  auxbasis = universal
  #PBE XC
  xName = 101
  cName = 130
  # FaT on
  frozen = False 
  nmax_FaT = 5
  super = False
  emb_error = True 
  method = 1
  ricore = 2000
  scfiterlimit = 30

Example 4 PbE+FaT for molecule-in-molecule embedding (supermolecular basis):

In order to perform an PbE+FaT calculation for molecular active and environment subsystems, the input file would look like the following

$EMBED
  nsystm = 2
  ptnIndx = 2 4
  #No KEDF
  KEfunc = electro
  basis = def2-SVP
  auxbasis = universal
  #PBE XC
  xName = 101
  cName = 130
  # FaT on
  frozen = False
  nmax_FaT = 5
  super = True
  emb_error = True 
  method = 3
  ricore = 2000
  scfiterlimit = 30

The error in PbE+FaT with respect to KS-DFT should be of the order of micro Hartrees.

Example 5 PbE+FaT for periodic-in-periodic embedding (supermolecular basis):

Consider a diamond cell with 3D periodicity and 8 atoms. The totalCoord looks like

$coord
  0.84253439282839      0.84253439282839      0.84253439282839      c
  5.89774074979874      5.89774074979874      5.89774074979874      c
  0.84253439282839      4.21267196414196      4.21267196414196      c
  4.21267196414196      0.84253439282839      4.21267196414196      c
  4.21267196414196      4.21267196414196      0.84253439282839      c
  5.89774074979874      2.52760317848518      2.52760317848518      c
  2.52760317848518      5.89774074979874      2.52760317848518      c
  2.52760317848518      2.52760317848518      5.89774074979874      c
$end

Let the first two atoms serve as active subssytem (periodic) and the remaining atoms as the environment subsystem (periodic). In order to perform a PbE+FaT calculation for periodic active and environment subsystems, the input file would look like the following

$EMBED
  nsystm = 2
  ptnIndx = 2 8
  KEfunc = electro
  basis = def2-SVP
  auxbasis = universal
  #PBE_XC
  xName = 101 
  cName = 130 
  frozen = False
  #frozen = True
  scratch = True
  super = True
  emb_error = True 
  method = 5
  periodicity = 3
  cell_params =   3.5668 3.5668 3.5668 90.0 90.0 90.0
  kpoints = nkpoints 10 10 10
  scfconv = 7
  denconv = 6
  ricore = 120000
  nmax_FaT = 5 
  mxitdiis = 10
  scfiterlimit = 100

Example 6 FDE for molecule-in-periodic embedding:

Consider an \(H_2\) molecule interacting with an \(H_{10}\) chain. The totalCoord file for this system is

$coord  natoms=    12
  6.29999918989539      0.69999990998838      3.75588736500000      h
  6.29999918989539     -0.69999990998838      3.75588736500000      h
  5.59999927990701      0.00000000000000      0.00000000000000      h
  6.99999909988376      0.00000000000000      0.00000000000000      h
  8.39999891986052      0.00000000000000      0.00000000000000      h
  0.00000000000000      0.00000000000000      0.00000000000000      h
  1.39999981997675      0.00000000000000      0.00000000000000      h
  2.79999963995350      0.00000000000000      0.00000000000000      h
  4.19999945993026      0.00000000000000      0.00000000000000      h
  9.79999873983727      0.00000000000000      0.00000000000000      h
 11.19999855981402      0.00000000000000      0.00000000000000      h
 12.59999837979077      0.00000000000000      0.00000000000000      h
$end

Let the \(H_2\) molecule as well as the two center atoms of the \(H_{10}\) chain be considered as the active subsystem (molecular) and the remaining as the environment (periodic). Molecule-in-periodic embedding using a KEDF-based embedding potential (as in FDE) can be performed by using an input file that looks like the following:

$EMBED
  nsystm = 2
  ptnIndx = 4 12
  KEfunc = 521
  basis = def2-TZVPP
  auxbasis = universal
  xName = 1
  cName = 12
  super = False 
  method = 2
  scfconv = 6
  denconv = 6
  ricore = 2000
  scfiterlimit = 40
  mxitdiis = 5
  cell_params = 7.408479999999999
  periodicity = 1
  kpoints = nkpoints 10

Please note, that the KEDF-based molecule-in-periodic embedding (method=2) works in a significantly different way than molecule-in-molecule embedding. The total energy of the embedded active subsystem does not include the contribution from the embedding potential. It should not be used to perform DFT-in-DFT embedding but rather WFT-in-DFT embedding. For more details of method 2 please refer to Sec. 2.6.2 of [213]

Example 7 WFT-in-DFT embedding:

After running any of the previous examples one gets a file called embdpot in the directory of the active subsystem (Cluster). Using this file one can perform WFT-in-DFT embedding for the active subsystem using the dscf and ricc2 modules. To do this the control file should contain the following text

$riperembdpot file=embdpot

Example 8 RT-TDDFET with frozen embedding potential:

Similar to WFT-in-DFT embedding, once the embedding potential for the molecular active subsystem is generated at the end of the riperembed run, one can use this embedding potential to perform RT-TDDFET with a frozen embedding potential by adding the following text to the control file and keeping the embdpot file in the working directory.

$fde
  activeSubsystem 1
  # super #Uncomment to use supermolecular basis
  method 4

Note, the user needs to also add the appropriate keywords and settings for the standard RT-TDDFT calculation. Finally, one needs to run riper.

Example 9 RT-TDDFET with updating embedding potential:

In the previous example, RT-TDDFET was performed with a frozen embedding potential. It is also possible to update the embedding potential partially, i.e., updating the terms corresponding to the active subsystem. Updating the embedding potential is supported only for method=1 (FDE) and method=3 (PbE). To do this one needs to firstly run the riperembed script with the appropriate input. Once this calculation is complete, one needs to go into the directory of the active subsystem (Cluster) and simply add the appropriate keywords and settings as for standard RT-TDDFT calculations and run riper.