OCT-MCTDH Documentation


Using optimal control theory together with MCTDH is realized through the Python script optcntrl that iteratively optimizes the electric field to achieve an optimization target.

The current implementation of optcntrl has been tested with Python versions 2.4.4 and 2.5.2 on Linux machines and assumes that the path to the Python executable is found in /usr/bin/env python. On some architectures this could be different. In this case the path can be changed in the first line of the file $MCTDH_DIR/bin/python/optcntrl.py.
NOTE: The OCT python scripts will not work with Python3. All other python scripts run on python3.

The script invokes the MCTDH main program for solving the time-dependent Schrödinger equation and a background process (efield) to evaluate the electric field during the propagation. Both programs are communicating through named pipes.

Starting an optimal control job

To start optcntrl on a platform where MCTDH is installed it type:
optcntrl  [-opt]  input_file
Possible options can be obtained with
optcntrl -h
Purpose    : Optimal control calculation by running MCTDH iteratively.

Usage      : optcntrl [ -mnd -c -D -n -r -R -w -x -deb -ver -h -? ] inputfile

  -mnd       Make name directory.
  -nd        Do not delete old output files in name directory before starting
             a new calculation.
  -D <dir>   Denotes the directory where files are written to,
             (name in .inp file ignored).
  -n <niter> Perform niter iterations (ignore value in input file)
  -r         Restart calculation. Perform additional iterations.
  -R         Restart calculation similar to -r, but start with a
             forward calculation (if crashed during forward).
  -w         Allow to overwrite existing data.
  -x <suf>   Use MCTDH/analyse executables with ending "suf",
             e.g., "-x .g95" to use binaries compiled with g95.
  -deb       Run in debug mode (more detailed logging)
  -ver       Print version info
  -h -?      Print this help text.

For restart-runs, the input file is parsed but no temporary input
files are created. The current iteration number is obtained from
file 'overlaps'.
To stop the calculation, create an empty file 'stop'.
I.e. type 'touch stop' while being in the working directory.
The argument needed by optcntrl is the name of an input file which is structured analogously to the input files of the MCTDH main Program. The input file is processed by the optcntrl script and temporary input files for the MCTDH main program are generated. Indeed, most parts of the input file remain unchanged and are used within temporary input files. Hence, it is possible to use most of the keywords and keyword placing as described in the MCTDH input documentation - exceptions are outlined below.

Input Documentation

Besides the Sections known from the MCTDH input two additional sections, the OCT-SECTION containing control parameters and the TARGET-SECTION specifying the control target, have to be given. An optional third section, the FILTER-SECTION, can be used to impose restrictions on the optimal field within the Krotov algorithm. Below only new sections and those that are modified compared to the original sections of the MCTDH main program are outlined. For other options refer to the MCTDH input documentation.


XXX Description
RUN  Propagation and output-times, sampling rate of the electrical field, etc
PRIMITIVE-BASIS Definition of primitive basis. See MCTDH input documentation
SPF-BASIS Definition of the single-particle function basis See MCTDH input documentation
OCT Control parameters, functional to be used, initial guess, etc.
FILTER Constraints on the optimal field
OPERATOR The operator to be used. See MCTDH input documentation
INIT_WF How to generate the initial wavefunction. See MCTDH input documentation
TARGET The target operator or target states to be used
INTEGRATOR The integrator to be used. See MCTDH input documentation


As within the MCTDH main program, the RUN-SECTION specifies the type of calculation to be performed, output to be written, etc. For an optimal control run a number of keywords are required:
Required keywords
propagation Indicate a calculation of type propagation
optcntrl(=S) Indicates an optimal control run. S can be either S=pc to run in a predictor corrector mode or S=simprop. In the latter case simultaneous forward and backward propagations of the initial and target states are performed. This may be useful if the wave function to be saved becomes very large. With S=simprop saving the wave function (i.e., psi=double) can be avoided. Note: if S is not set or S=pc then psi=double is required.
tfinal=R Sets the final time in fs where the control target has to be reached. Note that tinit is zero.
tout=R Sets the time step for the sampling of the electric field. Note: if tpsi is given, then it has to be equal to tout.

The name keyword is not required in the RUN-SECTION if the -D option is used. If name is given it specifies the work directory that contains all temporary files and output of the calculation.

All external files -- except the .op file containing the operator -- specified in the input files such as natpot files or files containing externally created initial and target states must be given either with absolute path or a path relative to the name directory.

As for the MCTDH program the flag overwrite allows overwriting of existing data in the name directory. This is the same as providing the -w option in the command line. Note, that existing files in the name directory are deleted if overwrite or -w is is given (except if the current job is restarted with the -r or -R option).

Exact calculations and distributed memory prallelization (MPI) are not implemented yet, i.e., the keywords exact and usempi must not be used.


The OCT-SECTION contains control parameters and parameters determining the initial guess field. Also the type of functional to be used is specified here.
algorithm = S S is the name of the OCT algorithm be used. S can be either default or Krotov. The default algorithm is the Rabitz OCT scheme.
In case of the Krotov scheme a number of additional keywords can be defined in the OCT section. See below.
guess = S S is the name of a function defined in the operator file specifying the guess field. If guess is not given, a zero guess field will be used.
penalty-function = S S is the name of a function defined in the operator file specifying a pulse envelope multiplied to the field. If penalty-function is not given, S will be set to unity.
reference = S S is the name of a function defined in the operator file specifying the reference field. If reference is not given, S will be set to zero. Note: the control functional is usually not monotonic convergent if S is some arbitrary time-dependent function. The reference keyword must not be used within the Krotov algorithm - here the reference-function is defined as the field from the previous iteration.
penalty = R R is the penalty-factor in atomic units. Default: R=1
iterations = I Sets the total number of iterations to be performed.
operator = S The operator S used as dipole operator.
write = S(,S1(,S2...)) Several output files can be created during the iterations. Possible options for S are
tddipmtx: Time-dependent matrix element of the dipole operator.
tdoverlap: Time-dependent overlap of initial and target states both at time t (requires target to be a wavefunction).
tdtargetpop: Time-dependent expectation value of the target.
rawfield: The new field before applying envelope, penalty-factor and reference-function.

Alternatively one can simply set write=all to write all possible output files.
save=S(,S1(,S2...)) Save output files as filename.i for all iterations i. save=Sj requires write=Sj. See write. In addition to the Sj given above, also save=field can be given (the field for the current iteration is written in any case.) Also save=all may be used to save all output.
tdoverlap Use time-dependent overlap instead of time-independent overlap to calculate the new field (requires target to be a wavefunction)
alignphs For multi-target runs optimize
1/N^2 |sum_i <Psi_i|Psi_tar,i>|^2 instead of 1/N sum_i |<Psi_i|Psi_tar,i>|^2, where N is the number of targets

The implementation of the Krotov OCT scheme offers a number of additional options as given in the table below.

filter (=S) Filtering of the optimal field. This option requires a FILTER-SECTION in the input file (see below). If S is not given, the filter defined in the FILTER-SECTION is interpreted as a stop-filter, i.e., as projection onto an unwanted subspace. This is also true for S=stop. If S=pass the filter defined in the FILTER-SECTION is interpreted as a pass-filter, i.e., as projection onto an wanted subspace. Additionally, if S=none, filtering is not applied and the FILTER-SECTION is ignored.

For details about filtering see Ref.

no-reffilter The reference-field remains unfiltered. See Ref.
save=S(,S1(,S2...)) As above but with the additional Keyword S=filter. If S=filter is given, temporary files used for filtering are saved for all iterations.


Within the Krotov scheme for obtaining the optimal field filtering of the field can be applied. The FILTER-SECTION is used to define a projector onto an (un-)wanted subspace of the complete search space of the field. The FILTER-SECTION is processed by the external program fldop so that in principle all operations available for program fldop can be used within the FILTER-SECTION. However, most of those operations are not projectors. At present, two keywords can be used to define a filter:

filter=S The field is Fourier-transformed and a function S defined in the operator file is multiplied to the transformed field in the frequency domain. Then the inverse transformation is applied. Note, that for a Projector S2=S has to be fullfiled and that the resulting field must be real, i.e. S has to be symmetric (S(ω)=S(-ω)).
symmetry=S Projection onto time-symmetric (S=even) or time-anti-symmetric (S=odd) fields.

NOTE: There is no check if the FILTER-SECTION defines a projector and if that projector is real-valued.


The target section can be used to either construct a target state, or to specify a target operator. If a target state is constructed, the TARGET-SECTION will be treated as an INIT_WF-SECTION. (see MCTDH input documentation)

If the expectation value of a target operator is to be optimized the only keyword given in this section is:

operator=S Optimize the expectation value of the operator S, where S is the label of an operator defined in the operator file. Note, that S has to be time-independent (OCT is not implemented for time-dependent target operators).

If operator is set, all other lines in this section are ignored.

Multi-Target optimizations

Multi target optimizations can be performed using the multi-packet algorithm. Within the SPF-BASIS-SECTION the keyword
has to be set, where n is the number of targets and initial states. Accordingly, within the INIT_WF-SECTION and the TARGET-SECTION multi-packet targets have to be constructed.

Note, that the multi-set algorithm (which is default) has to be used for multi-target optimizations. If the primitive basis includes electronic states, then also the states need to be propagated using the multi-set algorithm, i.e., both keywords

have to be set within the SPF-BASIS-SECTION.

Multiple targets can either be defined by defining a multi-packet wave function in the TARGET-SECTION or by defining a multi-packet operator. While the definition of wave functions is done streight forwardly, multi-target operators have to be defined in a single HAMILTONIAN-SECTION. The definition depends on wether or not an electronic basis is present in the PRIMITIVE-BASIS-SECTION. If no electronic basis is present one can define an additional mode called packets in the HAMILTONIAN-SECTION, e.g.:

  modes         |  X      | packets
  1             | rgtstep | S1&1
  1             | lftstep | S2&2
For packets, the same notation as for electronic states is used (since internally, packets are uncoupled electronic states). In the above example the first target may be a projection on the right side of the primitive grid of mode X while packet two is projected on the left side.

NOTE: In order to be able to address packets one has to set the usepack keyword in the HAMILTONIAN-SECTION.

If an electronic basis is present in the PRIMITIVE-BASIS-SECTION then the packet operator has to be constructed manually using the operator mode for electronic states. Suppose the system contains N states and P packets. Then the total number of internal electronic states is N*P. To address the single states one has to use the combined index of the N*P internal states. The nth state of the pth packet will have the comnined index N*(p-1)+n, i.e, the index of the electronic states defined in PRIMITIVE-BASIS-SECTION runs fastest followed by the package index.

The example below gives the target operator for a system containing two electronic states and two packets. Here, the target for packet one is the second electronic state while for packet two the target is the first electronic state.

  modes         |  el
  1             |  S2&2   # 2 = N*(p-1)+n with N=2, n=2, p=1
  1             |  S3&3   # 3 = N*(p-1)+n with N=2, n=1, p=2

Operator file

Definition of Operators

The operator file must contain at least two operators, the system operator - including the dipole operator multiplied with the electric field - and a separate operator containing the dipole operator alone. If the control target is to maximize the expectation value of an operator, e.g. the projector on an electronic state, this operator has to be defined as a separate operator as well.

Example input: system and dipole operator using a 2-D modified Henon-Heiles potential

# System Hamiltonian
  modes         |  X   | Y    | Time
1.0             |  KE  | 1    | 1
0.5             |  q^2 | 1    | 1
-lambda/3.0     |  q^3 | 1    | 1
lambda^2/16.0   |  q^4 | 1    | 1
1.0             |  1   | KE   | 1
0.5             |  1   | q^2  | 1
lambda^2/16.0   |  1   | q^4  | 1
lambda          |  q   | q^2  | 1
lambda^2/8.0    |  q^2 | q^2  | 1
 -mu            | 1    | q    | field
 mu/2.0         | q    | 1    | field

# Dipole operator as above but without electric field
    modes       |  X   | Y
 mu             | 1    | q
-mu/2.0         | q    | 1
The entry field in the Time column has to be specified in the LABELS-SECTION:
field = external1d{fieldf.dat}
The file fieldf.dat is generated by optcntrl and holds the sampled electrical field. Definition of analytic fields The scalar fields such as the guess-field and the pulse envelope have to be specified in the operator file. Similar to Hamiltonians analytic functions are defined with
  function terms
Example input:

 wr   = 800,cm-1
 tmax = 200,fs

 sin2=sin[1.0/tmax*PI/1.5 ,0.]^2

 sin2*sin1*0.01 +  sin2*0.01

All labels and parameters defined in the .op file can be used for the function definitions as well if they are one-dimensional analytic expressions (no derivatives etc.). In general, the terms in the function definitions follow the same rules as within a HAMILTON-SECTION, except that addition and substraction are allowed as well.

Note, that it is possible to use existing sampled data via the external1d statement in the LABELS-SECTION. In this case, the data files must contain data from -2*tout to tfinal+2*tout.

The PARAMETER-SECTION, LABELS-SECTION and the FUNCTION-SECTION are read by an external program genfld that generates an ASCII file with the sampled function with a sampling interval from 0 to tfinal and spacing tout. Note that tfinal is not known as a parameter to genfld so that it - if needed - has to be explicitly defined in the PARAMETER-SECTION. If used as a stand-alone tfinal has to be given as an additional parameter, e.g, using the -p option in the command line (see below).

Output files

Besides the usual output created by MCTDH several output files are written. The list below gives a short overview about the files and their structure. the suffix <i> either stands for a file written during the ith iteration or for 'dat' in which case the file was written during the last iteration.
  • Contains overlaps of initial and target states at the final time
  • Structure:
  • 	do iiter=0,maxiter
    	   iiter, (dble(ovlp(j)),imag(ovlp(j)),abs(ovlp(j))**2,j=1,ntarget)
    	end do
    iiter=0 is the overlap obtained with the guess-field.

  • Contains the control functional
  • Structure:
  • 	  do iiter=0,maxiter
    	     iiter, J1, J2, (J3,) J1 + J2 (- abs(J3))
    	  end do
    where J1 is the total target population and J2 is the contribution of the field (negative). If the Krotov algorithm with filtering is used then J3 contains the (signed) contribution resulting from the constraints. iiter=0 is the functional obtained with the guess-field.

fieldf.<i> / fieldb.<i>
  • Contains the field obtained during forward / backward iteration i
  • Structure:
  • 	  do i=-2,nfin+2
    	     i*tout, (fld(i,j), j=1, 3)
    	  end do
    where the 3rd and 4th columns contain the extrapolated fields. If optcntrl=pc was given in the RUN-SECTION then the 4th column contains the predicted field. Data is written in a.u.

  • Contains the target population obtained during a forward iteration i
  • Structure:
  • 	  do i=0,nfin
    	     i*tout,total_population, (ovtar(j), j=1, ntarget)
    	  end do
    where total_population is the total target population, i.e., J1 , and ovtar(j)=<Psi_j(t)|O_j|Psi_j(t)> of target operator O_j, with O_j being the projector onto the target state j. Data is written in a.u.

dipmtxf.<i> / dipmtxb.<i>
  • Contains the time-dependent matrix element mtx= <psi_tar(t)|dipole|psi_init(t)>
  • Structure:
  • 	  do i=0,nfin
    	  end do
    Data is written in a.u.

overlapf.<i> / overlapb.<i>
  • Contains the time-dependent overlap ovlp= <psi_tar(t)|psi_init(t)>
  • Structure:
  • 	  do i=0,nfin
    	  end do
    Data is written in a.u.

rawfieldf.<i> / rawfieldb.<i>
  • Contains the time-dependent 'raw' field rfld=-Im ovlp(t)*mtx(t)
  • Structure:
  • 	  do i=0,nfin
    	     i*tout, rfld
    	  end do
    Data is written in a.u.

Example Inputs



The program efield is used to calculate the new field during the iterations. It communicates with MCTDH through named pipes and is started in the background before MCTDH is started. This program is usually started by the optcntrl script.

 Purpose: Calculates the electric field for optimal control

 Usage: efieldver [-pc -b -A -T -O -F -S -R -p -i -e -d -t -m -o -r
 -ver -h -?]

 -pc      : Run in predictor-corrector mode.
 -b       : Compute the electric field during a backward propagation.
 -A       : Align phases of overlaps of multiple initial and target states
            i.e. optimize 1/N^2|sum_i <Psi_i|Targ,i>|^2
            instead of 1/N sum_i |<Psi_i|Targ,i>|^2
 -T oper  : Maximise matrix element <Psi|oper|Psi> of operator oper.
 -O oper  : Use oper as dipole operator to compute the field.
 -F pnlty : Parameter to penalize the field strength. Default: pnlty=1.
 -S       : Load penalty-function from file pnltyfnc.dat.
 -R       : Load reference-function from file reffnc.dat.
 -p       : Do not read old psi file. Instaed propagate 2 WFs simultaneously.
 -i iter  : Number of current iteration.
 -e num   : Use polynomial of order num (0,1,2) to extrapolate the field.
            Default: num=2.
 -d       : Use time-dependent overlap (it target is a WF).
 -t       : Write target population during a forward propagation.
 -m       : Write dipole matrix element as function of time.
 -o       : Write overlap with target as function of time.
 -r       : Write raw field -Im<tar|mu|psi> as function of time.
 -h or -? : Print this help text.
 -ver     : Print version information.

 The matrix element of the dipole operator (-O option) formed with the initial
 and target state is calculated and the E-field is evaluated.
 This routine must be started in the background before the invocation of MCTDH.
 It exchanges information with the running propagation via FIFOs.

 Option -O is mandatory. Options -T/-A and -p/-pc may not be used together.

 Note: This routine is usually started by the script optcntrl. Therefore all
 file names are hardcoded (see subroutine ef_open).


The program fldop is used to operate on existing one-dimensional fields. It can be used to modify existing fields within control algorithms, e.g. to impose constraints. At present only the filtering of the field is performed by fldop.

 Purpose:  Operate an operator P on a sampled 1D function.

 Usage  : fldop84 [-f -o -O -p -w -ver -h -?] inputfile

 -f ifile : File containing the initial sampled field.
 -o ofile : Name of the output file (rather then ifile1).
 -O file  : File containing function definitions (if needed)
 -C       : Operate 1-P on the field.
 -p S R   : Add real parameter with name S and value R to define the functions
            The value R may carry a unit. E.g.: -p en 5,eV
 -w       : Enable overwrite.
 -ver     : Print version information.
 -h, -?   : Print this help text.

 The option -f is mandatory.

 "inputfile" contains a list of operations enclosed in keywords
 "begin-fldop" and "end-fldop".
 All operations enclosed are applied multiplicative:
 the output the nth operation is used as input for the (n+1)th operation.
 All actions defined in the fldop - block therefore define the
 operator P.


The input file must contain a textblock of the form

  operation definitions
All other parts of the input file are ignored. The available operation definitions are listed in the table below. All operations have to be defined one per line and will be operated multiplicative, that is, the result of the previous operation is used as input for the next operation. Therefore the order in which the operations are given is relevant to the result.

multfnc=S The sampled field is multiplied with a function S defined in the operator file.
symmetry=S Projection onto time-symmetric (S=even) or anti-symmetric (S=odd) fields.
maxval=R The sampled field is scaled such that max(abs(field))=R
filter=R Apply Fourier-transform, multiply function S defined in operator file, apply inverse Fourier-transform. S is assumed to be defined in the frequency domain where the frequency ordering is -ωmax...0...+ωmax
norm=R Scale the sampled field such that the integral of field2=R
reverse reverse the ordering of the field, i.e., f(n)-> f(N-n+1), where N is the number of sampling points.
complex=S If S=real then the field is set to Re[field].
If S=imag then the field is set to Im[field].
If S=conjg then the field is set to its complex conjugated.
If S=abs then the field is set to its absolute value.
scale=R Multiply the sampled field with a real number R.
addfnc=S Add the function S defined in the operator file to the sampled field.


As mentioned above, the program genfld is used to generate one-dimensional functions used within the OCT-MCTDH script. It can, however, be used as a stand-alone. It requires an input file similar to the operator input files (.op). Here only the FUNCTION- and, if used, PARAMETER- and LABELS- sections are required. Other sections are ignored.

 Purpose: Generate a sampled 1D function (ASCII) from an analytic expression.

 Usage  : genfldver [-o -O -p -a -w -R -h -? -ver] tinit tfinal tout fncname

 -o file  : Name of the output file (default: fncname.dat)
 -O file  : File containing function definitions (.op file)
 -p S R   : Add real parameter with name S and value R to define the functions
            The value R may carry a unit. E.g.: -p en 5,eV
 -a num   : Write num data points before and after given sampling period
 -w       : Enable overwrite.
 -R       : Reverse, i.e., write f(T-t) instead of f(t)
 -ver     : Print version information.
 -h, -?   : Print this help text.

 The option -O is mandatory if fncname is not zero or unity.

 tinit and tfinal define the first and last sampling point,
 tout gives the time step, all in a.u.
 fncname is the label of a function defined in an operator file.
 fncname has the pre-defined labels zero and unity for which no operator
 file is needed. zero and unity can be overwritten in the .op file.
 Output is in a.u.


Purpose : Computation of the XFROG trace of a 1-dim function, using a
          rectangular gate function with Gaussian shoulders.

Usage   : xfrog [-o -n -rf -rt -fd -td -g -ver -h -?] infile

Options :
 -o FILE     : The output is written to file FILE rather than to infile.xfg
               The string FILE can be a relative or a full path-name.
 -n nfield   : Number of points to be read from infile (default: all)
 -rf nres    : Sampling of output in frequency domain (default: 200)
 -rt nres    : Sampling of output in time domain (default: 200)
 -fd min max : Frequency range of the output (default: 0 to nfield/T/10)
 -td min max : Time range of the output (default: 0 to T)
 -g w1 w2    : Width of the gatefunction in time domain (default: w1=0, w2=T/10)
               w1 width of the rectangular part,
               w2 width of the Gaussian shoulders.
 -ver        : Version info of this program.
 -h, -?      : Print this help text.

The argument infile is the file containing the electric field. The
gate function is rectangular pulse of unity height and width w1 with Gaussian
shoulders of the form exp(-t**2/(2*w2**2))
The user will  be prompted for missing arguments. All input and output is in au.


M. Schröder and A. Brown, "Generalized Filtering of laser Felds in optimal control theory: Application to symmetry Filtering of quantum gate operations", New Journal of Physics, Quantum control special issue, release date (tentative), October 2009.