Real-time TDDFT

Overview

Real-time time-dependent density functional theory (RT-TDDFT) is a DFT-based approach to electronic excited states based on integrating the time-dependent Kohn-Sham (TDKS) equations in time. The theoretical underpinnings, strengths, and limitations are similar to traditional linear-response (LR) TDDFT methods, but instead of a frequency domain solution to the TDKS equations, RT-TDDFT yields a full time-resolved, potentially non-linear solution. Real-time simulations can be used to compute not only spectroscopic properties (e.g., absorption spectra, polarizabilites, etc), but also the time and space-resolved electronic response to arbitrary external stimuli (e.g., electron charge dynamics after laser excitation). For theoretical and computational details, please refer to the following paper:

1. K. Lopata, N. Govind,”Modeling Fast Electron Dynamics with Real-Time Time-Dependent Density Functional Theory: Application to Small Molecules and Chromophores”, J. Chem. Theory Comput., 7, 1344 (2011) DOI:10.1021/ct200137z

This functionality is built on the Gaussian basis set DFT module, and will work for closed-shell (spin-restricted) and open-shell (spin unrestricted) calculations with essentially any combination of basis set and exchange-correlation functional in NWChem. The current implementation assumes frozen nuclei (Born-Oppenheimer approximation).

In a nutshell, running a RT-TDDFT calculation takes the following form:

1. Compute ground state density matrix with DFT module
2. Propagate density matrix using RT-TDDFT
3. Post-process resulting time-dependent observables (e.g., dipole moment)

Units

Unless specified otherwise, all inputs and outputs are in atomic units. Some useful conversions are:

Quantity	Conversion
Time	1 au = 0.02419 fs
Length	1 au = 0.5292 A
Energy	1 au = 27.2114 eV
Electric field	1 au = 514.2 V/nm
Dipole moment	1 au = 2.542 D

Syntax

The charge, geometry, basis set, and DFT options are all specified as normal, using their respective syntax. Real-time TDDFT parameters are supplied in the RT_TDDFT block (note, nothing is case-sensitive), with all possible options summarized below, and each discussed in detail afterwards.

RT_TDDFT 
  [TMAX <double default 1000>]  
  [DT <double default 0.1>]  
  [TAG <string default "<rt_tddft>: "] 
  [LOAD (scf || vectors <string>)]  
  [NCHECKS <integer default 10>] 
  [NPRINTS (* || <integer>)]  
  [NRESTARTS (* || <integer>)]  
  [TOLERANCES (zero <double default 1e-8> || series <double default 1e-10> || interpol <double default 1e-7>)]  
  [PROPAGATOR (euler || rk4 || magnus) default magnus]  
  [EXP (diag || pseries)]  
  [PROF]  
  [NOPROP]  
  [STATIC]  
  [PRINT (*, dipole, quadrupole, field, moocc, field, energy, cputime, charge, convergence, s2)] 
  [EXCITE <string geomname>` with <string fieldname>]  
  [FIELD] 
    ...  
  [END]  
  [VISUALIZATION]  
    ...  
  [END] 
  [LOAD RESTART]
END

TMAX: Simulation time

This option specifies the maximum time (in au) to run the simulation before stopping, which must be a positive real number. In practice, you can just stop the simulation early, so in most cases it is simplest to just set this to a large value to ensure you capture all the important dynamics (see the example of resonant ultraviolet excitation of water). For most valence excitations, for example, 1000 au is overkill so you might want to automatically stop at 500 au:

rt_tddft
  ...  
  tmax 500.0 
  ... 
end

DT: Time step

This specifies the electronic time step for time integration. A larger time step results in a faster simulation, but there are two issues which limit the size of the time step. First, the integration algorithm can become unstable and/or inaccurate for larger time steps, which depends on the integrator used. Very roughly speaking, the second order Magnus integrator should be reliable for time steps up to 0.2 au. Second, you must choose a time step small enough to capture the oscillations of interest, i.e., to resolve an excitation of energy ω, your time step needs to be smaller than π/ω, and typically a tenth of that for accuracy. For example, to capture high energy oscillations such as core-level excitations (e.g., ω = 50au) you might want a relative small time step:

rt_tddft
 ...
 dt 0.01
 ...
end

It is always good practice to check that your results are independent of time step.

TAG: Output label

This option sets a label for the output for convenient parsing (e.g., with “grep”). Every output line with time-dependent data will begin with this string (set to <rt_tddft>: by default). For example setting:

rt_tddft
 ...
 tag "nameofrun"
 ...
end

will result in outputs that look like:

...
nameofrun      2.20000  -7.589146713114E+001     # Etot
...

NCHECKS: Number of run-time check points

This option takes an integer number (default 10), or * which means every step, which sets the total of number of run-time checkpoints, where various sanity checks are made such as symmetries, idempotency, traces, etc. These checks are not terribly efficient (e.g., involve re-building the TD Fock matrix) so excessive checking will slow down execution.

rt_tddft
  ...
  nchecks 10
  ...
end

NPRINTS: Number of print points

This sets the number of print points, i.e., the total number of time-dependent observables (e.g., dipole, charge, energy) that are computed and printed. It either takes an integer number or * which means every time step (this is the default). Since there is no appreciable cost to computing and printing these quantities, there is usually no need to change this from *.

rt_tddft
  ...
  nprints *
  ...
end

NRESTARTS: Number of restart checkpoints

This sets the number of run-time check points where the time-dependent complex density matrix is saved to file, allowing the simulation to be restarted) from that point. By default this is set to 0. There is no significant computational cost to restart checkpointing, but of course there is some disk I/O cost (which may become somewhat significant for larger systems). For example, in the following example there will be 100 restart points, which corresponds to 1 backup every 100 time steps.

rt_tddft
  ...
  tmax 1000.0
  dt 0.1
  nrestarts 100
  ...
end

TOLERANCES: Controlling numerical tolerances

This option controls various numerical tolerances:

• zero: threshold for checks that quantities are zero, e.g., in symmetry checks (default 1e-8)
• series: numerical convergence for series, e.g., matrix exponentiation (default 1e-10)
• interpol: numerical convergence for interpolation, e.g., in Magnus propagator (default 1e-7)

Occasionally it is useful to loosen the interpolation tolerances if the Magnus interpolator requires an excessive amount of steps; usually this will not impact accuracy. For example, this sets the tolerances to their defaults with a looser interpolation tolerance:

rt_tddft
  ...
  tolerances zero 1d-8 series 1d-10 interpol 1e-5
end

PROPAGATOR: Selecting the integrator method

This selects the propagator (i.e., time integrator) method. Possible choices include euler for Euler integration (terrible, you should never use this), rk4 for 4th order Runge-Kutta, and magnus for 2nd order Magnus with self-consistent interpolation. In virtually all cases Magnus is superior in terms of stability. Euler or rk4 are perhaps only useful for debugging or simplicity (e.g., for code development).

rt_tddft
  ...
  propagator magnus
  ...
end

EXP: Selecting the matrix exponentiation method

This selects the method for exponentiation matrices. For now this can either be pseries for a contractive power series or diag for diagonalization. In general the power series (default) is faster.

rt_tddft
  ...
  exp diag
  ...
end

PROF: Run-time profiling

This turns on run-time profiling, which will display time spent in each component of the code (e.g., building components of the TD Fock matrix, properites, etc). This slows code down slightly and results in very verbose output.

rt_tddft
  ...
  prof
  ...
end

NOPROP: Skipping propagation

This options causes the RT-TDDFT module skip propagation, i.e., to initialize and finalize. For now this is largely useful for skipping to visualization post-processing without having to re-run a simulation.

rt_tddft
  ...
  noprop
  ...
end

STATIC: Force static Fock matrix

This option sets the static Fock matrix flag, meaning the time-dependent Fock matrix will not be recalculated at each time, but instead use the t=0 value. This will drastically increase the simulation speed since the bulk of the work is spent rebuilding the TD Fock matrix, but will give non-physical results. For example, using static to compute an absorption spectrum will result in excitations corresponding to the raw eigenvalue differences without electron-hole relaxation. This option has few uses besides dry-runs and debugging.

rt_tddft
  ...
  static
  ...
end

PRINT: Selecting time-dependent quantities to be printed

This sets the various time-dependent properties that are to be computed and printed at each print point. Note that for many of these options, the values are computed and printed for each geometry specified in the input deck, not only the active one (i.e., the one set using set geometry ... in the input deck). Possible choices are:

• dipole: Dipole moment
• quadrupole: Quadrupole moment
• field: External (applied) electric field
• moocc: Molecular orbital occupations
• energy: Components of system energy (e.g., core, XC, total, etc)
• cputime: CPU time taken in simulation so far (useful for checking scaling)
• charge: Electronic charge (computed from density matrix, not from the XC grid)
• convergence: Convergence information (e.g., from Magnus)
• s2: value (openshell only)
• *: Print all quantities

The defaults correspond to:

rt_tddft
  ...
  print dipole field energy convergence
  ...
end

FIELD: Sub-block for specifying external electric fields

This sub-block is used to specify external electric fields, each of which must be given a unique name. Numerous fields can be specified, but each will applied to the system only if an appropriate excitation rule is set. There are a few preset field types; others would have to be manually coded. Note the names are arbitrary, but chosen here to be descriptive:

field "kick"
  type delta      # E(t=0) = max; E(t>0) = 0
  polarization x  # = x, y, z
  max 0.0001      # maximum value of electric field
  spin total      # = alpha, beta, total (only valid for open-shell)
end

field "gpulse"
  type gaussian   # Gaussian enveloped quasi-monochromatic pulse: E(t) = max * exp( -(t-t0)^2 / 2s^2) * sin(w0*t + phi0)
  polarization x  # = x, y, z
  frequency 0.12  # frequency of laser pulse in au (e.g., 0.12 au = 3.27 eV)
  phase 0.0       # phase shift of laser pulse (in rad)
  center 200.0    # center of Gaussian envelope (in au time)
  width 50.0      # width of Gaussian pulse (in au time)
  max 0.0001      # maximum value of electric field
  spin total      # = alpha, beta, total (only valid for open-shell)
end

field "hpulse"
  type hann       # sin^2 (Hann) enveloped quasi-monochromatic pulse
  polarization x  # = x, y, z
  frequency 0.12  # frequency of laser pulse in au (e.g., 0.12 au = 3.27 eV)
  center 200.0    # center of Hann envelope (in au time)
  width 50.0      # width of Hann pulse (in au time)
  max 0.0001      # maximum value of electric field
  spin total      # = alpha, beta, total (only valid for open-shell)
end

field "resonant"
  type cw         # monochromatic continuous wave
  frequency 0.12  # frequency of laser pulse in au (e.g., 0.12 au = 3.27 eV)
  polarization x  # = x, y, z
  max 0.0001      # maximum value of electric field
  spin total      # = alpha, beta, total (only valid for open-shell)
end

EXCITE: Excitation rules

This sets the rules for applying external fields to the system. It takes the form excite <geom> with <field>, where <geom> is the name of a geometry fragment (defaults to “geometry” which is the default geometry name), and <field> is the name of a field structure. Assuming, for example, you have defined a field name kick this option takes the form (note that quotes are optional and shown for clarity):

rt_tddft
  ...
  excite "geometry" with "kick"
  ...
end

VISUALIZATION: Sub-block for controlling 3D visualization

This block is used to control visualization of the electronic charge density, which is typically used during a resonant excitation. This is a two stage process. During the propagation, a series of density matrices will be dumped to file (see options below). After propagation, if the dplot option is set, the code will read in options from a separate DPLOT block and convert the density matrix snapshots to a corresponding series of real-space charge density cube files.

visualization
  tstart 0.0        # start visualization at this time
  tend 100.0        # stop visualization at this time
  treference 0.0    # subtract density matrix at this time (useful for difference densities)
  dplot             # post-process density matrices into cube files after propagation
end

LOAD RESTART

This keyword needs to be added to restart a calculation. In the following example, the calculation will restart from the previous calculation and extend the run to the new tmax

rt_tddft
  ...
  tmax 10.0        # new end time
  load restart
end

MOCAP: Sub-block for molecular orbital complex absorbing potential

The K. Lopata and N. Govind, “Near and Above Ionization Electronic Excitations with Non-Hermitian Real-Time Time-Dependent Density Functional Theory”, Journal of Chemical Theory and Computation 9 (11), 4939-4946 (2013) DOI:10.1021/ct400569s

  mocap
    maxval 100.0     # clamp CAP at this value (in Ha)
    emin 0.5         # any MO with eigenvalue >= 0.5 Ha will have CAP applied to it
    prefac 1.0       # prefactor for exponential
    expconst 1.0     # exponential constant for CAP
    on               # turn on/off CAP
    nochecks         # disable checks for speed
    noprint          # don't print CAP value
  end

MAXVAL: Exponential Maximum Value

EMIN: Vacuum Energy Level

ON/OFF: Turn on/off CAP

The default value is off, i.e. no complex absorbing potential is applied.

Worked Examples

Absorption spectrum of water

Here, we compute the 6-31G/TD-PBE0 absorption spectrum of gas-phase water using RT-TDDFT. In the weak-field limit, these results are essentially identical to those obtained via traditional linear-response TDDFT. Although it involves significantly more work do use RT-TDDFT in this case, for very large systems with many roots a real-time approach becomes advantageous computationally and also does not suffer from algorithm stability issues.

To compute the absorption, we find the ground state of the system and then subject it to three delta-function excitations (x,y,z), which simultaneously excites all electronic modes of that polarization. The three resulting dipole moments are then Fourier transformed to give the frequency-dependent linear polarizability, and thus the absorption spectrum. The full input deck is RT_TDDFT_h2o_abs.nw and the corresponding output is RT_TDDFT_h2o_abs.nwo.gz.

title "Water TD-PBE0 absorption spectrum"
echo


start water

##                                                                                                                                                                           
## aug-cc-pvtz / pbe0 optimized                                                                                                                                              
##                                                                                                                                                                           
## Note: you are required to explicitly name the geometry                                                                                                                    
##                                                                                                                                                                           
geometry "system" units angstroms nocenter noautoz noautosym
  O     0.00000000    -0.00001441    -0.34824012
  H    -0.00000000     0.76001092    -0.93285191
  H     0.00000000    -0.75999650    -0.93290797
end

## Note: We need to explicitly set the "active" geometry even though there is only one geom.                                                                                 
set geometry "system"

## All DFT and basis parameters are inherited by the RT-TDDFT code                                                                                                           
basis
 * library 6-31G
end

dft
  xc pbe0
end

## Compute ground state of the system                                                                                                                                        
task dft energy

##                                                                                                                                                                           
## Now, we compute an x, y, and z kick simulation, which we give separate "tags" for post-processing.                                                                  
##                                                                                                                                                                           

unset rt_tddft:*
rt_tddft
  tmax 200.0
  dt 0.2

  tag "kick_x"

  field "kick"
    type delta
    polarization x
    max 0.0001
  end

  excite "system" with "kick"
 end
task dft rt_tddft

unset rt_tddft:*
rt_tddft
  tmax 200.0
  dt 0.2

  tag "kick_y"

  field "kick"
    type delta
    polarization y
    max 0.0001
  end

  excite "system" with "kick"
 end
task dft rt_tddft

unset rt_tddft:*
rt_tddft
  tmax 200.0
  dt 0.2

  tag "kick_z"

  field "kick"
    type delta
    polarization z
    max 0.0001
  end

  excite "system" with "kick"
 end
task dft rt_tddft

After running the simulation, we extract the x-dipole moment for the x-kick and similarly for the y and z-kicks (see “contrib/parsers” directory for this script or download here: RT_TDDFT_scripts.tgz ).

nw_rtparse.py -xdipole -px -tkick_x h2o_abs.nwo > x.dat
nw_rtparse.py -xdipole -py -tkick_y h2o_abs.nwo > y.dat
nw_rtparse.py -xdipole -pz -tkick_z h2o_abs.nwo > z.dat

Note, the syntax for extracting the x polarization for the x-kick, etc. Alternatively, we could grep and cut, or whatnot. This will give use the resulting time-dependent dipole moments:

Now, we need to take the Fourier transforms of these dipole moments to yield the the x,x element of the 3x3 linear polarizability tensor, and similarly for the y,y and z,z elements. Here I am using an FFT utility, although any discrete Fourier transform will do. To accelerate convergence of the FFT, I have damped the time signals by exp(-t /τ) which results in Lorentzians with FWHM of 2 / τ and have also “zero padded” the data with 50000 points. This is not critical for extracting frequencies, but creates “cleaner” spectra, although care must be taken to damp sufficiently if padding to avoid artifacts (see small ripples around 23 eV in plot below). After Fourier transform, I “paste” the three files together to easily plot the absorption, which is constructed from the trace of the polarizability matrix, i.e., the sum of the imaginary parts of the FFTs of the dipole moments.

S (ω) = (4πω)/(3cκ)Tr[Im α(ω)]

where c is the speed of light (137 in atomic units), κ is the kick electric field strength, and α(ω) is the linear polarizabilty tensor computed from the Fourier transforms of the time-dependent dipole moments. For example,

fft1d -d50 -z -p50000 < x.dat | rotate_fft > xw.dat
fft1d -d50 -z -p50000 < y.dat | rotate_fft > yw.dat
fft1d -d50 -z -p50000 < z.dat | rotate_fft > zw.dat
paste xw.dat yw.dat zw.day > s.dat

Here, you can just use your favorite Fourier transform utility or analysis software, but for convenience there is also a simple GNU Octave fft1d.m utility in the “contrib/parsers” directory of the trunk or download here: RT_TDDFT_scripts.tgz Note, options are hardcoded at the moment, so the switches above are not correct instead edit the file and run (also it reads file rather than redirect from stdin). Assuming the FFT output takes the form (w, Re, Im, Abs), to plot using gnuplot we would do:

gnuplot> plot "s.dat" u ($1*27.2114):($1*($3+$7+$11))

where we have scaled by 27.2114 to output in eV instead of atomic units, and we have not properly scaled to get the absolute oscillator strengths (thus our magnitudes are in “arbitrary units”). The real-time spectrum is shown below, along with the corresponding linear-response TDDFT excitations are shown in orange for comparison. Since we are in the weak field regime, the two are identical. Note the oscillator strengths are arbitrary and scaled, if not scaled the area under each RT-TDDFT curve should integrate to the linear response oscillator strength.

Resonant ultraviolet excitation of water

In this example we compute the time-dependent electron response to a quasi-monochromatic laser pulse tuned to a particular transition. We will use the results of the previous example (6-31G/PBE0 gas-phase water). First, we consider the absorption spectrum (computed previously) but plotted for the three polarizations (x,y,z) rather then as a sum. Say we are interested in the excitation near 10 eV. We can clearly see this is a z-polarized transition (green on curve). To selectively excite this we could use a continuous wave E-field, which has a delta-function, i.e., single frequency, bandwidth but since we are doing finite simulations we need a suitable envelope. The broader the envelope in time the narrower the excitation in frequency domain, but of course long simulations become costly so we need to put some thought into the choice of our envelope. In this case the peak of interest is spectrally isolated from other z-polarized peaks, so this is quite straightforward. The procedure is outlined below, and the corresponding frequency extent of the pulse is shown on the absorption figure in orange. Note that it only covers one excitation, i.e., the field selectively excites one mode. The full input deck is RT_TDDFT_h2o_resonant.nw and the output is RT_TDDFT_h2o_resonant.nwo.gz.

title "Water TD-PBE0 resonant excitation" 
echo
scratch_dir ./scratch
permanent_dir ./perm

start water

##
## aug-cc-pvtz / pbe0 optimized
##
## Note: you are required to explicitly name the geometry
##
geometry "system" units angstroms nocenter noautoz noautosym
  O     0.00000000    -0.00001441    -0.34824012
  H    -0.00000000     0.76001092    -0.93285191
  H     0.00000000    -0.75999650    -0.93290797
end

## Note: We need to explicitly set the "active" geometry even though there is only one geom.
set geometry "system" 

## All DFT and basis parameters are inherited by the RT-TDDFT code
basis
 * library 6-31G
end

dft
 xc pbe0
end

## Compute ground state of the system
task dft energy

##
## We excite the system with a quasi-monochromatic
## (Gaussian-enveloped) z-polarized E-field tuned to a transition at
## 10.25 eV.  The envelope takes the form:
##
## G(t) = exp(-(t-t0)^2 / 2s^2)
##
## The target excitation has an energy (frequency) of w = 0.3768 au
## and thus an oscillation period of T = 2 pi / w = 16.68 au
##
## Since we are doing a Gaussian envelope in time, we will get a
## Gaussian envelope in frequency (Gaussians are eigenfunctions of a
## Fourier transform), with width equal to the inverse of the width in
## time.  Say, we want a Gaussian in frequency with FWHM = 1 eV
## (recall FWHM = 2 sqrt (2ln2) s_freq) we want an s_freq = 0.42 eV =
## 0.0154 au, thus in time we need s_time = 1 / s_time = 64.8 au.
##
## Now we want the envelope to be effectively zero at t=0, say 1e-8
## (otherwise we get "windowing" effects).  Reordering G(t):
##
## t0 = t - sqrt(-2 s^2 ln G(t))
##
## That means our Gaussian needs to be centered at t0 = 393.3 au.
##
## The total simulation time will be 1000 au to leave lots of time to
## see oscillations after the field has passed.
##
rt_tddft
  tmax 1000.0
  dt 0.2

  field "driver"
    type gaussian
    polarization z
    frequency 0.3768  # = 10.25 eV
    center 393.3
    width 64.8
    max 0.0001
  end

  excite "system" with "driver"
 end
task dft rt_tddft

From the time-dependent dipole moment you can see the field driving the system into a superposition of the ground state and the one excited state, which manifests as monochromatic oscillations. After the field has passed the dipole oscillations continue forever as there is no damping in the system.

Charge transfer between a TCNE dimer

Here we compute the time-dependent charge oscillations between a TCNE (tetracyanoethylene) dimer separated by 3 Angstroms, where the top molecule starts neutral and the bottom one starts with a -1 charge. This somewhat non-physical starting point will lead to far-from-equilibrium dynamics as the charge violently sloshes between the two molecules, with the oscillation period a function of the molecular separation. The trick here is to use fragments by have multiple geometries in the input deck, where each fragment is converged separately, then assembled together without SCF to use as a starting point. We use a small but diffuse basis and a range-separated functional (CAM-B3LYP). The input deck is RT_TDDFT_tcne_dimer.nw and the full output is RT_TDDFT_tcne_dimer.nwo.gz.

title "Tetracyanoethylene dimer charge transfer"

echo
scratch_dir ./scratch
permanent_dir ./perm

start tcne
echo

##
## Each fragment optimized with cc-pvdz/B3LYP
##
geometry "bottom" units angstroms noautosym nocenter noautoz
 C    -1.77576486     0.66496556     0.00004199
 N    -2.94676621     0.71379797     0.00004388
 C    -0.36046718     0.62491168     0.00003506
 C     0.36049301    -0.62492429    -0.00004895
 C     1.77579907    -0.66504145    -0.00006082
 N     2.94680364    -0.71382258    -0.00006592
 C    -0.31262746    -1.87038951    -0.00011201
 N    -0.85519492    -2.90926164    -0.00016331
 C     0.31276207     1.87031662     0.00010870
 N     0.85498782     2.90938919     0.00016857
end

geometry "top" units angstroms noautosym nocenter noautoz
 C    -1.77576486     0.66496556     3.00004199
 N    -2.94676621     0.71379797     3.00004388
 C    -0.36046718     0.62491168     3.00003506
 C     0.36049301    -0.62492429     2.99995105
 C     1.77579907    -0.66504145     2.99993918
 N     2.94680364    -0.71382258     2.99993408
 C    -0.31262746    -1.87038951     2.99988799
 N    -0.85519492    -2.90926164     2.99983669
 C     0.31276207     1.87031662     3.00010870
 N     0.85498782     2.90938919     3.00016857
end


## dimer geometry is the union of bottom and top geometry
geometry "dimer" units angstroms noautosym nocenter noautoz
 C    -1.77576486     0.66496556     0.00004199
 N    -2.94676621     0.71379797     0.00004388
 C    -0.36046718     0.62491168     0.00003506
 C     0.36049301    -0.62492429    -0.00004895
 C     1.77579907    -0.66504145    -0.00006082
 N     2.94680364    -0.71382258    -0.00006592
 C    -0.31262746    -1.87038951    -0.00011201
 N    -0.85519492    -2.90926164    -0.00016331
 C     0.31276207     1.87031662     0.00010870
 N     0.85498782     2.90938919     0.00016857
#---
 C    -1.77576486     0.66496556     3.00004199
 N    -2.94676621     0.71379797     3.00004388
 C    -0.36046718     0.62491168     3.00003506
 C     0.36049301    -0.62492429     2.99995105
 C     1.77579907    -0.66504145     2.99993918
 N     2.94680364    -0.71382258     2.99993408
 C    -0.31262746    -1.87038951     2.99988799
 N    -0.85519492    -2.90926164     2.99983669
 C     0.31276207     1.87031662     3.00010870
 N     0.85498782     2.90938919     3.00016857
end


##
## C, N: 3-21++G
##
basis spherical
C    S
    172.2560000              0.0617669        
     25.9109000              0.3587940        
      5.5333500              0.7007130        
C    SP
      3.6649800             -0.3958970              0.2364600        
      0.7705450              1.2158400              0.8606190        
C    SP
      0.1958570              1.0000000              1.0000000        
C    SP
      0.0438000              1.0000000              1.0000000        
N    S
    242.7660000              0.0598657        
     36.4851000              0.3529550        
      7.8144900              0.7065130        
N    SP
      5.4252200             -0.4133010              0.2379720        
      1.1491500              1.2244200              0.8589530        
N    SP
      0.2832050              1.0000000              1.0000000        
N    SP
      0.0639000              1.0000000              1.0000000        
end


##
## Charge density fitting basis.
##
basis "cd basis"
C    S
      5.91553927E+02         0.31582020
      1.72117940E+02         0.87503863
      5.47992590E+01         2.30760524
C    S
      1.89590940E+01         1.0000000
C    S
      7.05993000E+00         1.0000000
C    S
      2.79484900E+00         1.0000000
C    S
      1.15863400E+00         1.0000000
C    S
      4.94324000E-01         1.0000000
C    S
      2.12969000E-01         1.0000000
C    P
      3.27847358E-01         1.0000000
C    P
      7.86833659E-01         1.0000000
C    P
      1.97101832E+00         1.0000000
C    D
      4.01330100E+00         1.0000000
C    D
      1.24750500E+00         1.0000000
C    D
      4.08148000E-01         1.0000000
C    F
      9.00000000E-01         1.0000000
N    S
      7.91076935E+02         0.41567506
      2.29450184E+02         1.14750694
      7.28869600E+01         3.01935767
N    S
      2.51815960E+01         1.0000000
N    S
      9.37169700E+00         1.0000000
N    S
      3.71065500E+00         1.0000000
N    S
      1.53946300E+00         1.0000000
N    S
      6.57553000E-01         1.0000000
N    S
      2.83654000E-01         1.0000000
N    P
      4.70739194E-01         1.0000000
N    P
      1.12977407E+00         1.0000000
N    P
      2.83008403E+00         1.0000000
N    D
      5.83298650E+00         1.0000000
N    D
      1.73268650E+00         1.0000000
N    D
      5.45242500E-01         1.0000000
N    F
      1.82648000E+00         1.0000000
end


##
## Universal DFT parameters. Note, we are doing open-shell even for
## the neutral fragment so the movecs have the correct size. 
##
## We are using the CAM-B3LYP functional (no need to use "direct"
## since we are doing CD fitting).
##
dft
  xc xcamb88 1.00 lyp 0.81 vwn_5 0.19 hfexch 1.00
  cam 0.33 cam_alpha 0.19 cam_beta 0.46
  odft
  convergence density 1d-9
  grid fine
  maxiter 1000
end


##
## Converge bottom fragment with extra electron and top fragment as
## neutral.
##
charge -1
set geometry "bottom"
dft
  mult 2
  vectors input atomic output "bottom.movecs"
end
task dft energy

charge 0
set geometry "top"
dft
  mult 1
  vectors input atomic output "top.movecs"
end
task dft energy


##
## Assemble the two fragments but don't do SCF--this keeps the system
## in a far-from-equilibrium state from which we will watch the
## dynamics.
##
charge -1
set geometry "dimer"
dft
  mult 2
  vectors input fragment "bottom.movecs" "top.movecs" output "dimer.movecs"
  noscf
end
task dft energy


##
## Now do RT-TDDFT from this crazy state without any electric fields.
##
rt_tddft
  tmax 500.0
  dt 0.2
  load vectors "dimer.movecs"

  print dipole field energy s2 charge
end
task dft rt_tddft

The time-dependent charge shows that the excess electron starts on the “bottom” molecule (i.e., a total electronic charge of -65), then swings to completely occupy the “top” molecule then oscillates back and forth. The frequency of this oscillation is dependent on the separation, with larger separations leading to lower frequencies. It is important to note, however, this starting point is highly non-physical, specifically converging the two fragments together and “gluing” them together introduces an indeterminate amount of energy to the system, but this simulation shows how charge dynamics simulations can be done.

MO CAP example

## aug-cc-pvtz/PBE0 optimized
geometry "system" units angstroms noautosym noautoz nocenter
  O     0.00000043     0.11188833     0.00000000
  H     0.76000350    -0.47275229     0.00000000
  H    -0.76000393    -0.47275063     0.00000000
end

set geometry "system"

basis spherical
  * library 6-31G*
end

dft
  xc pbe0
  convergence density 1d-9
end
task dft

rt_tddft
  dt 0.2
  tmax 250.0

  print dipole field energy charge

  mocap
    expconst 1.0     # exponential constant for CAP
    emin 0.5         # any MO with eigenvalue >= 0.5 Ha will have CAP applied to it
    prefac 1.0       # prefactor for exponential
    maxval 100.0     # clamp CAP at this value (in Ha)
    on               # turn on CAP
    nochecks         # disable checks for speed
    noprint          # don't print CAP value
  end

  field "kick"
    type delta
    max 0.0001
    polarization z
  end

  excite "system" with "kick"
end
task dft rt_tddft

After running the parser script on the output files from the input above using either with the on or off keywords, the following plot can be produced from the data file obtained with the command nw_rtparse.py -xdipole -pz -t"<rt_tddft>"