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
timedependent Schrödinger equation and a background process
(efield) to evaluate the electric field during the propagation.
Both programs are communicating through named pipes.
Possible options can be obtained with
optcntrl [opt] input_file
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 restartruns, 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.
XXX  Description 

RUN  Propagation and outputtimes, sampling rate of the electrical field, etc 
PRIMITIVEBASIS  Definition of primitive basis. See MCTDH input documentation 
SPFBASIS  Definition of the singleparticle 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 


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 RUNSECTION 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.


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. 
penaltyfunction = S  S is the name of a function defined in the operator file specifying a pulse envelope multiplied to the field. If penaltyfunction 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 timedependent function. The reference keyword must not be used within the Krotov algorithm  here the referencefunction is defined as the field from the previous iteration. 
penalty = R  R is the penaltyfactor 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: Timedependent matrix element of the dipole operator. tdoverlap: Timedependent overlap of initial and target states both at time t (requires target to be a wavefunction). tdtargetpop: Timedependent expectation value of the target. rawfield: The new field before applying envelope, penaltyfactor and referencefunction. 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 timedependent overlap instead of timeindependent overlap to calculate the new field (requires target to be a wavefunction) 
alignphs  For multitarget runs optimize 1/N^2 sum_i <Psi_iPsi_tar,i>^2 instead of 1/N sum_i <Psi_iPsi_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 FILTERSECTION in
the input file (see below). If S is not given, the filter defined
in the FILTERSECTION is interpreted as a stopfilter, i.e., as
projection onto an unwanted subspace. This is also true for
S=stop. If S=pass the filter defined in the
FILTERSECTION is interpreted as a passfilter, i.e., as projection
onto an wanted subspace. Additionally, if S=none, filtering
is not applied and the FILTERSECTION is ignored. For details about filtering see Ref. 
noreffilter  The referencefield 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. 


filter=S  The field is Fouriertransformed 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 S^{2}=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 timesymmetric (S=even) or timeantisymmetric (S=odd) fields. 
NOTE: There is no check if the FILTERSECTION defines a projector and if that projector is realvalued.
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 timeindependent (OCT is not implemented for timedependent target operators). 
If operator is set, all other lines in this section are ignored.
has to be set, where n is the number of targets and initial states. Accordingly, within the INIT_WFSECTION and the TARGETSECTION multipacket targets have to be constructed.
packets=n
Note, that the multiset algorithm (which is default) has to be used for multitarget optimizations. If the primitive basis includes electronic states, then also the states need to be propagated using the multiset algorithm, i.e., both keywords
have to be set within the SPFBASISSECTION.
packets=n multiset
Multiple targets can either be defined by defining a multipacket wave function in the TARGETSECTION or by defining a multipacket operator. While the definition of wave functions is done streight forwardly, multitarget operators have to be defined in a single HAMILTONIANSECTION. The definition depends on wether or not an electronic basis is present in the PRIMITIVEBASISSECTION. If no electronic basis is present one can define an additional mode called packets in the HAMILTONIANSECTION, e.g.:
HAMILTONIANSECTION_proj usepack  modes  X  packets  1  rgtstep  S1&1 1  lftstep  S2&2 endhamiltoniansectionFor 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 HAMILTONIANSECTION.
If an electronic basis is present in the PRIMITIVEBASISSECTION 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*(p1)+n, i.e, the index of the electronic states defined in PRIMITIVEBASISSECTION 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.
HAMILTONIANSECTION_swap usepack  modes  el  1  S2&2 # 2 = N*(p1)+n with N=2, n=2, p=1 1  S3&3 # 3 = N*(p1)+n with N=2, n=1, p=2 endhamiltoniansection
Example input: system and dipole operator using a 2D modified HenonHeiles potential
# System Hamiltonian HAMILTONIANSECTION  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  endhamiltoniansection # Dipole operator as above but without electric field HAMILTONIANSECTION_dipole  modes  X  Y  mu  1  q mu/2.0  q  1  endhamiltoniansectionThe entry field in the Time column has to be specified in the LABELSSECTION:
LABELSSECTION field = external1d{fieldf.dat} ENDLABELSSECTIONThe file fieldf.dat is generated by optcntrl and holds the sampled electrical field.
Example input:
FUNCTIONSECTION_name function terms ENDFUNCTIONSECTION
PARAMETERSECTION wr = 800,cm1 tmax = 200,fs endparametersection LABELSSECTION sin2=sin[1.0/tmax*PI/1.5 ,0.]^2 sin1=sin[wr,0.0] endlabelssection FUNCTIONSECTION_guess sin2*sin1*0.01 + sin2*0.01 endfunctionsection FUNCTIONSECTION_envelope sin2 endfunctionsectionAll labels and parameters defined in the .op file can be used for the function definitions as well if they are onedimensional analytic expressions (no derivatives etc.). In general, the terms in the function definitions follow the same rules as within a HAMILTONSECTION, 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 LABELSSECTION. In this case, the data files must contain data from 2*tout to tfinal+2*tout.
The PARAMETERSECTION, LABELSSECTION and the FUNCTIONSECTION 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 PARAMETERSECTION. If used as a standalone tfinal has to be given as an additional parameter, e.g, using the p option in the command line (see below).
do iiter=0,maxiter iiter, (dble(ovlp(j)),imag(ovlp(j)),abs(ovlp(j))**2,j=1,ntarget) end doiiter=0 is the overlap obtained with the guessfield.
do iiter=0,maxiter iiter, J1, J2, (J3,) J1 + J2 ( abs(J3)) end dowhere 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 guessfield.
do i=2,nfin+2 i*tout, (fld(i,j), j=1, 3) end dowhere the 3rd and 4th columns contain the extrapolated fields. If optcntrl=pc was given in the RUNSECTION then the 4th column contains the predicted field. Data is written in a.u.
do i=0,nfin i*tout,total_population, (ovtar(j), j=1, ntarget) end dowhere total_population is the total target population, i.e., J1 , and ovtar(j)=<Psi_j(t)O_jPsi_j(t)> of target operator O_j, with O_j being the projector onto the target state j. Data is written in a.u.
do i=0,nfin i*tout,(dble(mtx(j)),imag(mtx(j)),j=1,ntarget) end doData is written in a.u.
do i=0,nfin i*tout,(dble(ovlp(j)),imag(ovlp(j)),j=1,ntarget) end doData is written in a.u.
do i=0,nfin i*tout, rfld end doData is written in a.u.
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 ?] Options: pc : Run in predictorcorrector 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^2sum_i <Psi_iTarg,i>^2 instead of 1/N sum_i <Psi_iTarg,i>^2 T oper : Maximise matrix element <PsioperPsi> 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 penaltyfunction from file pnltyfnc.dat. R : Load referencefunction 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 timedependent 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<tarmupsi> 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 Efield 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 onedimensional 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 Options: 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 1P 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 "beginfldop" and "endfldop". 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
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.
beginfldop operation definitions endfldop


multfnc=S  The sampled field is multiplied with a function S defined in the operator file. 
symmetry=S  Projection onto timesymmetric (S=even) or antisymmetric (S=odd) fields. 
maxval=R  The sampled field is scaled such that max(abs(field))=R 
filter=R  Apply Fouriertransform, multiply function S defined in operator file, apply inverse Fouriertransform. 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 field^{2}=R 
reverse  reverse the ordering of the field, i.e., f(n)> f(Nn+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 onedimensional functions used within the OCTMCTDH script. It can, however, be used as a standalone. 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 Options: 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(Tt) 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 predefined 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 1dim 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 pathname. 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.