Time-dependent propagation

OK, it is about time to do TDDFT. Let us perform a calculation of the time evolution of the density of the methane molecule.

Ground-state

The reason to do first a ground-state DFT calculation is that this will be the initial state for the real-time calculation.

Run the input file below to obtain a proper ground-state (stored in the restart directory) as from the total energy convergence tutorial. (Now we are using a more accurate bond length than before.)


CalculationMode = gs
UnitsOutput = eV_Angstrom

Radius = 3.5*angstrom
Spacing = 0.18*angstrom

CH = 1.097*angstrom
%Coordinates
  "C" |           0 |          0 |           0
  "H" |  CH/sqrt(3) | CH/sqrt(3) |  CH/sqrt(3)
  "H" | -CH/sqrt(3) |-CH/sqrt(3) |  CH/sqrt(3)
  "H" |  CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
  "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%

Time-dependent

Input

Now that we have a starting point for our time-dependent calculation, we modify the input file so that it looks like this:


CalculationMode = td
UnitsOutput = eV_Angstrom

Radius = 3.5*angstrom
Spacing = 0.18*angstrom

CH = 1.097*angstrom
%Coordinates
  "C" |           0 |          0 |           0
  "H" |  CH/sqrt(3) | CH/sqrt(3) |  CH/sqrt(3)
  "H" | -CH/sqrt(3) |-CH/sqrt(3) |  CH/sqrt(3)
  "H" |  CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
  "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%

Tf  = 0.1/eV
dt = 0.002/eV

TDPropagator = aetrs
TDMaxSteps = Tf/dt
TDTimeStep = dt


PoissonSolver = isf

Here we have switched on a new run mode, CalculationMode = td instead of gs, to do a time-dependent run. We have also added three new input variables:

Note that, for convenience, we have previously defined a couple of variables, Tf and dt. We have made use of one of the possible propagators, aetrs. The manual explains about the possible options; in practice this choice is usually good except for very long propagation where the etrs propagator can be more stable.

Output

Now run Octopus. This should take a few seconds, depending on the speed of your machine. The most relevant chunks of the standard output are



Input: [IonsConstantVelocity = no]
Input: [Thermostat = none]
Input: [MoveIons = no]
Input: [TDIonicTimeScale = 1.000]
Input: [TDTimeStep = 0.2000E-02 hbar/eV]
Input: [TDPropagationTime = 0.1000 hbar/eV]
Input: [TDMaxSteps = 50]
Input: [TDDynamics = ehrenfest]
Input: [TDScissor = 0.000]
Input: [TDPropagator = aetrs]
Input: [TDExponentialMethod = taylor]

Input: [IonsConstantVelocity = no]
Input: [Thermostat = none]
Input: [MoveIons = no]
Input: [TDIonicTimeScale = 1.000]
Input: [TDTimeStep = 0.2000E-02 hbar/eV]
Input: [TDPropagationTime = 0.1000 hbar/eV]
Input: [TDMaxSteps = 50]
Input: [TDDynamics = ehrenfest]
Input: [TDScissor = 0.000]
Input: [TDPropagator = aetrs]
Input: [TDExponentialMethod = taylor]


********************* Time-Dependent Simulation **********************
  Iter           Time        Energy   SC Steps    Elapsed Time
**********************************************************************

      1       0.002000   -218.785065         1         0.049
      2       0.004000   -218.785065         1         0.047
      3       0.006000   -218.785065         1         0.047
...

     48       0.096000   -218.785065         1         0.047
     49       0.098000   -218.785065         1         0.047
     50       0.100000   -218.785065         1         0.049

It is worthwhile to comment on a few things:

The time step

A key parameter is, of course, the time step, TDTimeStep. Before making long calculations, it is worthwhile spending some time choosing the largest time-step possible, to reduce the number of steps needed. This time-step depends crucially on the system under consideration, the spacing, on the applied perturbation, and on the algorithm chosen to approximate the evolution operator.

In this example, try to change the time-step and to rerun the time-evolution. Make sure you are using TDMaxSteps of at least 100, as with shorter runs an instability might not appear yet. You will see that for time-steps larger than 0.0024 the propagation gets unstable and the total energy of the system is no longer conserved. Very often it diverges rapidly or even becomes NaN.

Also, there is another input variable that we did not set explicitly, relying on its default value, TDExponentialMethod. Since most propagators rely on algorithms to calculate the action of the exponential of the Hamiltonian, one can specify which algorithm can be used for this purpose.

You may want to learn about the possible options that may be taken by TDExponentialMethod, and TDPropagator – take a look at the manual. You can now try some exercises:

Laser fields

Now we will add a time-dependent external perturbation (a laser field) to the molecular Hamiltonian. For brevity, we will omit the beginning of the file, as this is left unchanged. The relevant part of the input file, with the modifications and additions is:


Tf  = 1/eV
dt = 0.002/eV

TDPropagator = aetrs
TDMaxSteps = Tf/dt
TDTimeStep = dt

FromScratch = yes

amplitude = 1*eV/angstrom
omega = 18.0*eV
tau0 = 0.5/eV
t0 = tau0

%TDExternalFields
  electric_field | 1 | 0 | 0 | omega | "envelope_cos"
%

%TDFunctions
  "envelope_cos" | tdf_cosinoidal | amplitude | tau0 | t0
%

%TDOutput
 laser
 multipoles
%
Laser field used in the tutorial
Laser field used in the tutorial

The most important variables here are the TDExternalFields block and the associated TDFunctions block. You should carefully read the manual page dedicated to these variables: the particular laser pulse that we have employed is the one whose envelope function is a cosine.

Now you are ready to set up a run with a laser field. Be careful to set a total time of propagation able to accommodate the laser shot, or even more if you want to see what happens afterwards. You may also want to consult the meaning of the variable TDOutput.

A couple of important comments:


  1. A. Castro, M.A.L. Marques, and A. Rubio, Propagators for the time-dependent Kohn–Sham equations, J. Chem. Phys. 121 3425-3433 (2004); ↩︎