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.
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 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.
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 |
|
|
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.
|
|
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. |
|
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.
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.
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.
packets=n
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.
packets=n multi-set
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.:
HAMILTONIAN-SECTION_proj usepack --------------------------------------- modes | X | packets --------------------------------------- 1 | rgtstep | S1&1 1 | lftstep | S2&2 end-hamiltonian-sectionFor 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.
HAMILTONIAN-SECTION_swap usepack --------------------------------------- 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 end-hamiltonian-section
Example input: system and dipole operator using a 2-D modified Henon-Heiles potential
# System Hamiltonian HAMILTONIAN-SECTION --------------------------------------- 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 --------------------------------------- end-hamiltonian-section # Dipole operator as above but without electric field HAMILTONIAN-SECTION_dipole --------------------------------------- modes | X | Y --------------------------------------- mu | 1 | q -mu/2.0 | q | 1 --------------------------------------- end-hamiltonian-sectionThe entry field in the Time column has to be specified in the LABELS-SECTION:
LABELS-SECTION field = external1d{fieldf.dat} END-LABELS-SECTIONThe file fieldf.dat is generated by optcntrl and holds the sampled electrical field.
Example input:
FUNCTION-SECTION_name function terms END-FUNCTION-SECTION
PARAMETER-SECTION wr = 800,cm-1 tmax = 200,fs end-parameter-section LABELS-SECTION sin2=sin[1.0/tmax*PI/1.5 ,0.]^2 sin1=sin[wr,0.0] end-labels-section FUNCTION-SECTION_guess sin2*sin1*0.01 + sin2*0.01 end-function-section FUNCTION-SECTION_envelope sin2 end-function-sectionAll 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).
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 guess-field.
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 guess-field.
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 RUN-SECTION 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_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.
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 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 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 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
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.
begin-fldop operation definitions end-fldop
|
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 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(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.