How To Calculate and Write Out Expectation Values

Observables in the config file

Obtaining time-dependent expectation values for arbitrary observables while propagating a system is possible by defining observables in the config file. For example, when propagating a wave packet on a single potential surface, we may define the following `observables section <../config_file.html#observables>`__ in the config file.

observables: is_real = T, outfile = exp_values.dat, time_unit = ns
* type = r,   column_label='<r>', square='<r^2>', exp_unit = nm
* type = mom, column_label='<p>', square='<p^2>'

After the propagation (i.e. a call to `prop <../reference/prop.html#prop>`__), the file exp_values.dat will have been created and contain 5 columns: the time, and the expectation values of the position operator, the square of the position operator, the momentum operator, and the square of the momentum operator.

The expectation values of the squares may be used to later calculate the variance of the position or momentum operator. Whether or not these values are calculated and written out is determined entirely by the presence of the square parameter (which gives the column label). If square is not given, the expectation value of the square is not calculated, and the outfile will only contain a single column for that observable.

It is possible to write different expectation values to different files by specifying a different outfile for each observable. It is recommended to group all observables whose expectation value go into the same file into one section in the config file.

A large variety of observables may be defined in the config file, through the type parameter. This includes r and mom as in the above example, as well as ham, and kin to calculate energy expectation values (total and kinetic). Moreover, any grid-based operator may be specified, using the same parameters as in the `ham section <../config_file.html#ham>`__. For example, the following section would request the expectation values of the total, kinetic, and potential energy of a simple harmonic oscillator:

observables: is_real = T, outfile = exp_energies.dat
* type = ham,   column_label='<H>'
* type = kin,   column_label='<T>'
* type = op,    column_label='<V>', &
  op_type = pot, op_form = hamos, op_surf = 1, w_0 = 10

The most general definition of an observable is through type=matrix, which allows to define an arbitrary operator over the entire Hilbert space, read as a sparse matrix from file (cf. type=matrix in the `ham section <../config_file.html#ham>`__ of the config file).

Lastly, type=norm and type=pop are also available. The former calculates the norm of the wave function, the latter must be paired with a specification of exp_surf and calculates the population on the given energy surface.

The exp_surf parameter can be used more generally to calculate expectation values on a single energy surface, e.g.

observables: is_real = T, outfile = exp_values.dat, time_unit = ns, &
& type = r, exp_unit = nm
* column_label='<r> (tot)'
* column_label='<r> (sf1)', square='<r^2> (sf1)', exp_surf = 1
* column_label='<r> (sf2)', square='<r^2> (sf2)', exp_surf = 2

In this case, before calculating the expectation value, the wave function is projected onto the given surface and then renormalized.

Observables can be non-Hermitian, resulting in complex expectation values. This is possible for type=ham if the Hamiltonian is non-Hermitian, or with type=matrix, which allows to define completely arbitrary matrices. In this case, is_real = F must be specified. The output file will then contain two columns for the observable, one for the real part of the expectation value, and one for the imaginary part. If square is also defined, the number of columns for that observable will increase to a total of four.

The information about which expectation values should be calculated is determined from the config file data by the `init_prop routine <../reference/prop.html#init_prop>`__. In case any observable section in the config file uses the label parameter, only the items that match gen%ham%label in init_prop are taken into consideration. The information is stored as part of the prop_work variable.

The calculation and output of the expectation values are handled by the `prop routine <../reference/prop.html#prop>`__. Every time prop is called, the outfile for any of the observables is overwritten. While the propagation is running, data is continuously written to each outfile. Thus, running tail -f {outfile} may be a useful way to monitor the progress of a long-running propagation. On completion of the prop routine, the expectation values for the entire time grid covered by the propagation will be contained in outfile. For a backwards propagation the time (first column in outfile) will also run backwards.

Accessing expectation values programatically

In addition to writing expectation values to file, the `prop routine <../reference/prop.html#prop>`__ also stores all values in memory, in the complex matrix prop_work%exp%values. The number of rows of the matrix is equal to the number of time grid points. The first columns of the matrix contains the time grid values. The subsequent columns contain the expectation values for all observables, in the order in which they appear in the config file. Any observable that defines a square label has an additional column for the expectation value of the square of the observable. A label for all the columns may be obtained from the array prop_work%exp%column_labels.

Note that since prop_work%exp%values is complex, any observable with is_real=F does not result in two columns in the matrix. When processing columns that are known to contain the expectation values of Hermitian operators, the values may be converted to real on demand. All values are stored in internal units.

Expectation values will be included in prop_work%exp%values even if the observable does not define an outfile. In this case, the expectation value will not be written to file, but will be available in memory. It is even possible to not define an outfile for any of the observables, in order to not generate any output files, but only calculate expectation values internally.

The content of prop_work%exp%values is overwritten whenever the prop routine is called. Any values from an earlier call are lost.

Using prop_info_hook

While in most cases, the output of expectation values can entirely be specified in the config file, it can be useful to calculate some custom quantities using the prop_info_hook parameter of the `prop routine <../reference/prop.html#prop>`__. The parameter specifies a routine that is called after each propagation step. Inside the routine, expectation values could be calculated using the routines defined in the `exp module <../reference/exp.html>`__ (use qdyn_exp_mod). However, the most common usage of the prop_info_hook is to write out the wave function after each propagation step. E.g., using an external routine

subroutine prop_write_psi(psi, grid, ham, ti, work, para)

use qdyn
use qdyn_exp_mod
implicit none
type(state_t),     intent(inout) :: psi
type(grid_t),      intent(in)    :: grid
type(ham_t),       intent(inout) :: ham
integer,           intent(in)    :: ti
type(prop_work_t), intent(inout) :: work
type(para_t),      intent(in)    :: para

real(idp) :: t, r
type(state_t), save :: state_temp
character(len=datline_l) :: comment

t = para%tgrid%t_start + ti * para%tgrid%dt

call mold_state(state_temp, state)
call exp_r(r, state, state_temp, grid)

comment='# t = '//trim(nr2str(t))//'; <r> = '//trim(nr2str(r))
if (ti == 0) then
    call write_psi(psi, grid, 1, 1, 'psi1.dat',  &
    &              append=.false., comment=comment)
else
    call write_psi(psi, grid, 1, 1, 'psi1.dat',  &
    &              append=.true., comment=comment)
end if

end subroutine prop_write_psi

the propagation may be called as

call prop(psi, grid, gen, work, para, pulses, prop_info_hook=prop_write_psi)