**CalculationMode**
*Section*: Calculation Modes
*Type*: integer
*Default*: gs

Decides what kind of calculation is to be performed.
*Options*:

**gs**: Calculation of the ground state.**unocc**: Calculation of unoccupied/virtual KS states. Can also be used for a non-self-consistent calculation of states at arbitrary k-points, if`density.obf`from`gs`is provided in the`restart/gs`directory.**td**: Time-dependent calculation (experimental for periodic systems).**go**: Optimization of the geometry.**opt_control**: Optimal control.**em_resp**: Calculation of the electromagnetic response: electric polarizabilities and hyperpolarizabilities and magnetic susceptibilities (experimental for periodic systems).**casida**: Excitations via Casida linear-response TDDFT; for finite systems only.**vdw**: Calculate van der Waals coefficients.**vib_modes**: Calculation of the vibrational modes.**one_shot**: Obsolete. Use`gs`with`MaximumIter = 0`instead.**kdotp**: Calculation of effective masses by \(\vec{k} \cdot \vec{p}\) perturbation theory (experimental).**dummy**: This calculation mode does nothing. Useful for debugging, testing and benchmarking.**invert_ks**: Invert the Kohn-Sham equations (experimental).**recipe**: Prints out a tasty recipe.

**GOCenter**
*Section*: Calculation Modes::Geometry Optimization
*Type*: logical
*Default*: no

(Experimental) If set to yes, Octopus centers the geometry at
every optimization step. It also reduces the degrees of
freedom of the optimization by using the translational
invariance.

**GOConstrains**
*Section*: Calculation Modes::Geometry Optimization
*Type*: block

If `XYZGOConstrains`, `PDBConstrains`, and `XSFGOConstrains`
are not present, `Octopus` will try to fetch the geometry optimization
contrains from this block. If this block is not present, `Octopus`
will not set any constrains. The format of this block can be
illustrated by this example:

`%GOConstrains
'C' | 1 | 0 | 0
'O' | 1 | 0 | 0
%`

Coordinates with a constrain value of 0 will be optimized, while coordinates with a constrain different from zero will be kept fixed. So, in this example the x coordinates of both atoms will remain fixed and the distance between the two atoms along the x axis will be constant.

Note: It is important for the constrains to maintain the ordering in which the atoms were defined in the coordinates specifications. Moreover, constrains impose fixed absolute coordinates, therefore constrains are not compatible with GOCenter = yes

**GOFireIntegrator**
*Section*: Calculation Modes::Geometry Optimization
*Type*: integer
*Default*: verlet

The Fire algorithm (`GOMethod = fire`) uses a molecular dynamics
integrator to compute new geometries and velocities.
Currently, two integrator schemes can be selected
*Options*:

**euler**: The Euler method.**verlet**: The Velocity Verlet algorithm.

**GOFireMass**
*Section*: Calculation Modes::Geometry Optimization
*Type*: float
*Default*: 1.0 amu

The Fire algorithm (`GOMethod = fire`) assumes that all degrees of freedom
are comparable. All the velocities should be on the same
scale, which for heteronuclear systems can be roughly
achieved by setting all the atom masses equal, to the value
specified by this variable.
By default the mass of a proton is selected (1 amu).
However, a selection of `GOFireMass = 0.01` can, in manys systems,
speed up the geometry optimization procedure.
If `GOFireMass` <= 0, the masses of each
species will be used.

**GOLineTol**
*Section*: Calculation Modes::Geometry Optimization
*Type*: float
*Default*: 0.1

Tolerance for line-minimization. Applies only to GSL methods
that use the forces.
WARNING: in some weird units.

**GOMaxIter**
*Section*: Calculation Modes::Geometry Optimization
*Type*: integer
*Default*: 200

Even if the convergence criterion is not satisfied, the minimization will stop
after this number of iterations.

**GOMethod**
*Section*: Calculation Modes::Geometry Optimization
*Type*: integer
*Default*: fire

Method by which the minimization is performed. For more information see the
GSL documentation.
*Options*:

**steep_native**: (Experimental) Non-gsl implementation of steepest descent.**steep**: Simple steepest descent.**cg_fr**: Fletcher-Reeves conjugate-gradient algorithm. The conjugate-gradient algorithm proceeds as a succession of line minimizations. The sequence of search directions is used to build up an approximation to the curvature of the function in the neighborhood of the minimum.**cg_pr**: Polak-Ribiere conjugate-gradient algorithm.**cg_bfgs**: Vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) conjugate-gradient algorithm. It is a quasi-Newton method which builds up an approximation to the second derivatives of the function*f*using the difference between successive gradient vectors. By combining the first and second derivatives, the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.**cg_bfgs2**: The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher,*Practical Methods of Optimization*, Algorithms 2.6.2 and 2.6.4.**simplex**: This is experimental, and in fact,**not**recommended unless you just want to fool around. It is the Nead-Melder simplex algorithm, as implemented in the GNU Scientific Library (GSL). It does not make use of the gradients (*i.e.*, the forces) which makes it less efficient than other schemes. It is included here for completeness, since it is free.**fire**: The FIRE algorithm. See also`GOFireMass`and`GOFireIntegrator`. Ref: E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, and P. Gumbsch,*Phys. Rev. Lett.***97**, 170201 (2006).

**GOMinimumMove**
*Section*: Calculation Modes::Geometry Optimization
*Type*: float

Convergence criterion, for stopping the minimization. In
units of length; minimization is stopped when the coordinates
of all species change less than `GOMinimumMove`, or the
`GOTolerance` criterion is satisfied.
If `GOMinimumMove < 0`, this criterion is ignored.
Default is -1, except 0.001 b with `GOMethod = simplex`.
Note that if you use `GOMethod = simplex`,
then you must supply a non-zero `GOMinimumMove`.

**GOObjective**
*Section*: Calculation Modes::Geometry Optimization
*Type*: integer
*Default*: minimize_energy

This rather esoteric option allows one to choose which
objective function to minimize during a geometry
minimization. The use of this variable may lead to
inconsistencies, so please make sure you know what you are
doing.
*Options*:

**minimize_energy**: Use the total energy as objective function.**minimize_forces**: Use \(\sqrt{\sum_i \left| f_i \right|^2}\) as objective function. Note that in this case one still uses the forces as the gradient of the objective function. This is, of course, inconsistent, and may lead to very strange behavior.

**GOStep**
*Section*: Calculation Modes::Geometry Optimization
*Type*: float

Initial step for the geometry optimizer. The default is 0.5.
WARNING: in some weird units.
For the FIRE minimizer, default value is 0.1 fs,
and corresponds to the initial time-step for the MD.

**GOTolerance**
*Section*: Calculation Modes::Geometry Optimization
*Type*: float
*Default*: 0.001 H/b (0.051 eV/A)

Convergence criterion, for stopping the minimization. In
units of force; minimization is stopped when all forces on
ions are smaller than this criterion, or the
`GOMinimumMove` is satisfied. If `GOTolerance < 0`,
this criterion is ignored.

**PDBGOConstrains**
*Section*: Calculation Modes::Geometry Optimization
*Type*: string

Like `XYZGOConstrains` but in PDB format, as in `PDBCoordinates`.

**XSFGOConstrains**
*Section*: Calculation Modes::Geometry Optimization
*Type*: string

Like `XYZGOConstrains` but in XCrySDen format, as in `XSFCoordinates`.

**XYZGOConstrains**
*Section*: Calculation Modes::Geometry Optimization
*Type*: string

`Octopus` will try to read the coordinate-dependent constrains from the XYZ file
specified by the variable `XYZGOConstrains`.
Note: It is important for the contrains to maintain the ordering
in which the atoms were defined in the coordinates specifications.
Moreover, constrains impose fixed absolute coordinates, therefore
constrains are not compatible with GOCenter = yes

**InvertKSConvAbsDens**
*Section*: Calculation Modes::Invert KS
*Type*: float
*Default*: 1e-5

Absolute difference between the calculated and the target density in the KS
inversion. Has to be larger than the convergence of the density in the SCF run.

**InvertKSGodbyMu**
*Section*: Calculation Modes::Invert KS
*Type*: float
*Default*: 1.0

prefactor for iterative KS inversion convergence scheme from Godby based on van Leeuwen scheme

**InvertKSGodbyPower**
*Section*: Calculation Modes::Invert KS
*Type*: float
*Default*: 0.05

power to which density is elevated for iterative KS inversion convergence
scheme from Godby based on van Leeuwen scheme

**InvertKSMaxIter**
*Section*: Calculation Modes::Invert KS
*Type*: integer
*Default*: 200

Selects how many iterations of inversion will be done in the iterative scheme

**InvertKSStellaAlpha**
*Section*: Calculation Modes::Invert KS
*Type*: float
*Default*: 0.05

prefactor term in iterative scheme from L Stella

**InvertKSStellaBeta**
*Section*: Calculation Modes::Invert KS
*Type*: float
*Default*: 1.0

residual term in Stella iterative scheme to avoid 0 denominators

**InvertKSTargetDensity**
*Section*: Calculation Modes::Invert KS
*Type*: string
*Default*: `target_density.dat`

Name of the file that contains the density used as the target in the
inversion of the KS equations.

**InvertKSVerbosity**
*Section*: Calculation Modes::Invert KS
*Type*: integer
*Default*: 0

Selects what is output during the calculation of the KS potential.
*Options*:

**0**: Only outputs the converged density and KS potential.**1**: Same as 0 but outputs the maximum difference to the target density in each iteration in addition.**2**: Same as 1 but outputs the density and the KS potential in each iteration in addition.

**InvertKSmethod**
*Section*: Calculation Modes::Invert KS
*Type*: integer
*Default*: iterative

Selects whether the exact two-particle method or the iterative scheme
is used to invert the density to get the KS potential.
*Options*:

**two_particle**: Exact two-particle scheme.**iterative**: Iterative scheme for \(v_s\).**iter_stella**: Iterative scheme for \(v_s\) using Stella and Verstraete method.**iter_godby**: Iterative scheme for \(v_s\) using power method from Rex Godby.

**KSInversionAsymptotics**
*Section*: Calculation Modes::Invert KS
*Type*: integer
*Default*: xc_asymptotics_none

Asymptotic correction applied to \(v_{xc}\).
*Options*:

**xc_asymptotics_none**: Do not apply any correction in the asymptotic region.**xc_asymptotics_sc**: Applies the soft-Coulomb decay of \(-1/\sqrt{r^2+1}\) to \(v_{xc}\) in the asymptotic region.

**KSInversionLevel**
*Section*: Calculation Modes::Invert KS
*Type*: integer
*Default*: ks_inversion_adiabatic

At what level `Octopus` shall handle the KS inversion.
*Options*:

**ks_inversion_none**: Do not compute KS inversion.**ks_inversion_adiabatic**: Compute exact adiabatic \(v_{xc}\).

**OCTCheckGradient**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 0.0

When doing QOCT with the conjugate-gradient optimization scheme, the gradient is
computed thanks to a forward-backwards propagation. For debugging purposes, this
gradient can be compared with the value obtained "numerically" (*i.e.* by doing
successive forward propagations with control fields separated by small finite
differences).

In order to activate this feature, set `OCTCheckGradient` to some non-zero value,
which will be the finite difference used to numerically compute the gradient.

**OCTClassicalTarget**
*Section*: Calculation Modes::Optimal Control
*Type*: block

If `OCTTargetOperator = oct_tg_classical`, the you must supply this block.
It should contain a string (e.g. "(q[1,1]-q[1,2])*p[2,1]") with a mathematical
expression in terms of two arrays, q and p, that represent the position and momenta
of the classical variables. The first index runs through the various classical particles,
and the second index runs through the spatial dimensions.

In principle, the block only contains one entry (string). However, if the expression is very
long, you can split it into various lines (one column each) that will be concatenated.

The QOCT algorithm will attempt to maximize this expression, at the end of the propagation.

**OCTControlFunctionOmegaMax**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: -1.0

The Fourier series that can be used to represent the control functions must be truncated;
the truncation is given by a cut-off frequency which is determined by this variable.

**OCTControlFunctionRepresentation**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: control_fourier_series_h

If `OCTControlRepresentation = control_function_parametrized`, one must
specify the kind of parameters that determine the control function.
If `OCTControlRepresentation = control_function_real_time`, then this variable
is ignored, and the control function is handled directly in real time.
*Options*:

**control_fourier_series_h**: The control function is expanded as a full Fourier series (although it must, of course, be a real function). Then, the total fluence is fixed, and a transformation to hyperspherical coordinates is done; the parameters to optimize are the hyperspherical angles.**control_zero_fourier_series_h**: The control function is expanded as a Fourier series, but assuming (1) that the zero frequency component is zero, and (2) the control function, integrated in time, adds up to zero (this essentially means that the sum of all the cosine coefficients is zero). Then, the total fluence is fixed, and a transformation to hyperspherical coordinates is done; the parameters to optimize are the hyperspherical angles.**control_fourier_series**: The control function is expanded as a full Fourier series (although it must, of course, be a real function). The control parameters are the coefficients of this basis-set expansion.**control_zero_fourier_series**: The control function is expanded as a full Fourier series (although it must, of course, be a real function). The control parameters are the coefficients of this basis-set expansion. The difference with the option`control_fourier_series`is that (1) that the zero-frequency component is zero, and (2) the control function, integrated in time, adds up to zero (this essentially means that the sum of all the cosine coefficients is zero).**control_rt**: (experimental)

**OCTControlFunctionType**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: controlfunction_mode_epsilon

The control function may fully determine the time-dependent form of the
external field, or only the envelope function of this external field, or its phase.
Or, we may have two different control functions, one of them providing the phase
and the other one, the envelope.

Note that, if `OCTControlRepresentation = control_function_real_time`, then the control
function must **always** determine the full external field (THIS NEEDS TO BE FIXED).
*Options*:

**controlfunction_mode_epsilon**: In this case, the control function determines the full control function: namely, if we are considering the electric field of a laser, the time-dependent electric field.**controlfunction_mode_f**: The optimization process attempts to find the best possible envelope. The full control field is this envelope times a cosine function with a "carrier" frequency. This carrier frequency is given by the carrier frequency of the`TDExternalFields`in the`inp`file.

**OCTCurrentFunctional**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: oct_no_curr

(Experimental) The variable `OCTCurrentFunctional` describes which kind of
current target functional \(J1_c[j]\) is to be used.
*Options*:

**oct_no_curr**: No current functional is used, no current calculated.**oct_curr_square**: Calculates the square of current \(j\): \(J1_c[j] = {\tt OCTCurrentWeight} \int{\left| j(r) \right|^2 dr}\). For`OCTCurrentWeight`< 0, the current will be minimized (useful in combination with target density in order to obtain stable final target density), while for`OCTCurrentWeight`> 0, it will be maximized (useful in combination with a target density in order to obtain a high-velocity impact, for instance). It is a static target, to be reached at total time.**oct_max_curr_ring**: Maximizes the current of a quantum ring in one direction. The functional maximizes the \(z\) projection of the outer product between the position \(\vec{r}\) and the current \(\vec{j}\): \(J1[j] = {\tt OCTCurrentWeight} \int{(\vec{r} \times \vec{j}) \cdot \hat{z} dr}\). For`OCTCurrentWeight`> 0, the current flows in counter-clockwise direction, while for`OCTCurrentWeight`< 0, the current is clockwise.**oct_curr_square_td**: The time-dependent version of`oct_curr_square`. In fact, calculates the square of current in time interval [`OCTStartTimeCurrTg`, total time =`TDMaximumIter`*`TDTimeStep`]. Set`TDPropagator`=`crank_nicolson`.

**OCTCurrentWeight**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 0.0

In the case of simultaneous optimization of density \(n\) and current \(j\), one can tune the importance
of the current functional \(J1_c[j]\), as the respective functionals might not provide results on the
same scale of magnitude. \(J1[n,j]= J1_d[n]+ {\tt OCTCurrentWeight}\ J1_c[j]\). Be aware that its
sign is crucial for the chosen `OCTCurrentFunctional` as explained there.

**OCTDelta**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 0.0

If `OCTScheme = oct_mt03`, then you can supply the "eta" and "delta" parameters
described in [Y. Maday and G. Turinici, *J. Chem. Phys.* **118**, 8191 (2003)], using the
`OCTEta` and `OCTDelta` variables.

**OCTDirectStep**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 0.25

If you choose `OCTScheme = oct_direct` or `OCTScheme = oct_nlopt_bobyqa`,
the algorithms necessitate an initial "step" to perform the direct search for the
optimal value. The precise meaning of this "step" differs.

**OCTDoubleCheck**
*Section*: Calculation Modes::Optimal Control
*Type*: logical
*Default*: true

In order to make sure that the optimized field indeed does its job, the code
may run a normal propagation after the optimization using the optimized field.

**OCTDumpIntermediate**
*Section*: Calculation Modes::Optimal Control
*Type*: logical
*Default*: true

Writes to disk the laser pulse data during the OCT algorithm at intermediate steps.
These are files called `opt_control/laser.xxxx`, where `xxxx` is the iteration number.

**OCTEps**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 1.0e-6

Define the convergence threshold. It computes the difference between the "input"
field in the iterative procedure, and the "output" field. If this difference is
less than `OCTEps` the iteration is stopped. This difference is defined as:

\(
D[\varepsilon^{in},\varepsilon^{out}] = \int_0^T dt \left| \varepsilon^{in}(t)-\varepsilon^{out}(t)\right|^2
\)

(If there are several control fields, this difference is defined as the sum over
all the individual differences.)

Whenever this condition is satisfied, it means that we have reached a solution point
of the QOCT equations, *i.e.* a critical point of the QOCT functional (not
necessarily a maximum, and not necessarily the global maximum).

**OCTEta**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 1.0

If `OCTScheme = oct_mt03`, then you can supply the "eta" and "delta" parameters
described in [Y. Maday and G. Turinici, *J. Chem. Phys.* **118**, 8191 (2003)], using the
`OCTEta` and `OCTDelta` variables.

**OCTExcludedStates**
*Section*: Calculation Modes::Optimal Control
*Type*: string

If the target is the exclusion of several targets, ("OCTTargetOperator = oct_exclude_states")
then you must declare which states are to be excluded, by setting the OCTExcludedStates variable.
It must be a string in "list" format: "1-8", or "2,3,4-9", for example. Be careful to include
in this list only states that have been calculated in a previous "gs" or "unocc" calculation,
or otherwise the error will be silently ignored.

**OCTFilter**
*Section*: Calculation Modes::Optimal Control
*Type*: block

The block `OCTFilter` describes the type and shape of the filter function
that are applied to the optimized laser field in each iteration.
The filter forces the laser field to obtain the given form in frequency space.
Each line of the block describes a filter; this way you can actually have more
than one filter function (*e.g.* a filter in time and two in frequency space).
The filters are applied in the given order, *i.e.*, first the filter specified
by the first line is applied, then second line.
The syntax of each line is, then:

`%OCTFilter
domain | function
%`

Possible arguments for domain are:

(i)

(ii)

Example:

time | "exp(-80*( w + 0.1567 )^2 ) + exp(-80*( w - 0.1567 )^2 )"

%

Be careful that also the negative-frequency component is filtered since the resulting field has to be real-valued.

**frequency_filter**: The filter is applied in the frequency domain.

**OCTFixFluenceTo**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 0.0

The algorithm tries to obtain the specified fluence for the laser field.
This works only in conjunction with either the WG05 or the straight iteration scheme.

If this variable is not present in the input file, by default the code will not
attempt a fixed-fluence QOCT run. The same holds if the value given to this
variable is exactly zero.

If this variable is given a negative value, then the target fluence will be that of
the initial laser pulse given as guess in the input file. Note, however, that
first the code applies the envelope provided by the `OCTLaserEnvelope` input
option, and afterwards it calculates the fluence.

**OCTFixInitialFluence**
*Section*: Calculation Modes::Optimal Control
*Type*: logical
*Default*: yes

By default, when asking for a fixed-fluence optimization (`OCTFixFluenceTo = whatever`),
the initial laser guess provided in the input file is scaled to match this
fluence. However, you can force the program to use that initial laser as the initial
guess, no matter the fluence, by setting `OCTFixInitialFluence = no`.

**OCTHarmonicWeight**
*Section*: Calculation Modes::Optimal Control
*Type*: string
*Default*: "1"

(Experimental) If `OCTTargetOperator = oct_tg_plateau`, then the function to optimize is the integral of the
harmonic spectrum \(H(\omega)\), weighted with a function \(f(\omega)\)
that is defined as a string here. For example, if
you set `OCTHarmonicWeight = "step(w-1)"`, the function to optimize is
the integral of \(step(\omega-1)*H(\omega)\), *i.e.*
\(\int_1^{\infty} H \left( \omega \right) d\omega\).
In practice, it is better if you also set an upper limit, *e.g.*
\(f(\omega) = step(\omega-1) step(2-\omega)\).

**OCTInitialState**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: oct_is_groundstate

Describes the initial state of the quantum system.
Possible arguments are:
*Options*:

**oct_is_groundstate**: Start in the ground state.**oct_is_excited**: Currently not in use.**oct_is_gstransformation**: Start in a transformation of the ground-state orbitals, as defined in the block`OCTInitialTransformStates`.**oct_is_userdefined**: Start in a user-defined state.

**OCTInitialTransformStates**
*Section*: Calculation Modes::Optimal Control
*Type*: block

If `OCTInitialState = oct_is_gstransformation`, you must specify an
`OCTInitialTransformStates` block, in order to specify which linear
combination of the states present in `restart/gs` is used to
create the initial state.

The syntax is the same as the `TransformStates` block.

**OCTInitialUserdefined**
*Section*: Calculation Modes::Optimal Control
*Type*: block

Define an initial state. Syntax follows the one of the `UserDefinedStates` block.
Example:

`%OCTInitialUserdefined
1 | 1 | 1 | "exp(-r^2)*exp(-i*0.2*x)"
%`

**OCTLaserEnvelope**
*Section*: Calculation Modes::Optimal Control
*Type*: block

Often a pre-defined time-dependent envelope on the control function is desired.
This can be achieved by making the penalty factor time-dependent.
Here, you may specify the required time-dependent envelope.

It is possible to choose different envelopes for different control functions.
There should be one line for each control function. Each line should
have only one element: a string with the name of a time-dependent function,
that should be correspondingly defined in a `TDFunctions` block.

**OCTLocalTarget**
*Section*: Calculation Modes::Optimal Control
*Type*: string

If `OCTTargetOperator = oct_tg_local`, then one must supply a function
that defines the target. This should be done by defining it through a string, using
the variable `OCTLocalTarget`.

**OCTMaxIter**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: 10

The maximum number of iterations.
Typical values range from 10-100.

**OCTMomentumDerivatives**
*Section*: Calculation Modes::Optimal Control
*Type*: block

This block should contain the derivatives of the expression given in
`OCTClassicalTarget` with respect to the p array components.
Each line corresponds to a different classical particle, whereas the
columns correspond to each spatial dimension: the (i,j) block component
corresponds with the derivative wrt p[i,j].

**OCTNumberCheckPoints**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: 0

During an OCT propagation, the code may write the wavefunctions at some time steps (the
"check points"). When the inverse backward or forward propagation
is performed in a following step, the wavefunction should reverse its path
(almost) exactly. This can be checked to make sure that it is the case -- otherwise
one should try reducing the time-step, or altering in some other way the
variables that control the propagation.

If the backward (or forward) propagation is not retracing the steps of the previous
forward (or backward) propagation, the code will write a warning.

**OCTOptimizeHarmonicSpectrum**
*Section*: Calculation Modes::Optimal Control
*Type*: block
*Default*: no

(Experimental)
If `OCTTargetOperator = oct_tg_hhg`, the target is the harmonic emission spectrum.
In that case, you must supply an `OCTOptimizeHarmonicSpectrum` block in the `inp`
file. The target is given, in general, by:

\(J_1 = \int_0^\infty d\omega \alpha(\omega) H(\omega)\),

where \(H(\omega)\) is the harmonic spectrum generated by the system, and
\(\alpha(\omega)\) is some function that determines what exactly we want
to optimize. The role of the `OCTOptimizeHarmonicSpectrum` block is to determine
this \(\alpha(\omega)\) function. Currently, this function is defined as:

\(\alpha(\omega) = \sum_{L=1}^{M} \frac{\alpha_L}{a_L} \sqcap( (\omega - L\omega_0)/a_L )\),

where \(\omega_0\) is the carrier frequency. \(M\) is
the number of columns in the `OCTOptimizeHarmonicSpectrum` block. The values of *L* will be listed
in the first row of this block; \(\alpha_L\) in the second row, and \(a_L\) in
the third.

Example:

`%OCTOptimizeHarmonicSpectrum
7 | 9 | 11
-1 | 1 | -1
0.01 | 0.01 | 0.01
%`

**OCTPenalty**
*Section*: Calculation Modes::Optimal Control
*Type*: float
*Default*: 1.0

The variable specifies the value of the penalty factor for the
integrated field strength (fluence). Large value = small fluence.
A transient shape can be specified using the block `OCTLaserEnvelope`.
In this case `OCTPenalty` is multiplied with time-dependent function.
The value depends on the coupling between the states. A good start might be a
value from 0.1 (strong fields) to 10 (weak fields).

Note that if there are several control functions, one can specify this
variable as a one-line code, each column being the penalty factor for each
of the control functions. Make sure that the number of columns is equal to the
number of control functions. If it is not a block, all control functions will
have the same penalty factor.

All penalty factors must be positive.

**OCTPositionDerivatives**
*Section*: Calculation Modes::Optimal Control
*Type*: block

This block should contain the derivatives of the expression given in
`OCTClassicalTarget` with respect to the q array components.
Each line corresponds to a different classical particle, whereas the
columns correspond to each spatial dimension: the (i,j) block component
corresponds with the derivative wrt q[i,j].

**OCTRandomInitialGuess**
*Section*: Calculation Modes::Optimal Control
*Type*: logical
*Default*: false

The initial field to start the optimization search is usually given in the `inp` file,
through a `TDExternalFields` block. However, you can start from a random guess if you
set this variable to true.

Note, however, that this is only valid for the "direct" optimization schemes; moreover
you still need to provide a `TDExternalFields` block.

**OCTScheme**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: OPTION__OCTSCHEME__ZBR98

Optimal Control Theory can be performed with `Octopus` with a variety of different
algorithms. Not all of them can be used with any choice of target or control function
representation. For example, some algorithms cannot be used if
`OCTControlRepresentation = control_function_real_time`
(`OCTScheme` > `oct_straight_iteration`), and others cannot be used
if `OCTControlRepresentation = control_function_parametrized`
(`OCTScheme` < `oct_straight_iteration`).
*Options*:

**oct_nlopt_bobyqa**: The BOBYQA algorithm, as implemented in the NLOPT library -- therefore, octopus has to be compiled with it in order to be able to use this option.**oct_nlopt_lbfgs**: The local BFGS, as implemented in the NLOPT library -- therefore, octopus has to be compiled with it in order to be able to use this option.**oct_zbr98**: Backward-Forward-Backward scheme described in*JCP***108**, 1953 (1998). Only possible if target operator is a projection operator. Provides fast, stable and monotonic convergence.**oct_zr98**: Forward-Backward-Forward scheme described in*JCP***109**, 385 (1998). Works for projection and more general target operators also. The convergence is stable but slower than ZBR98. Note that local operators show an extremely slow convergence. It ensures monotonic convergence.**oct_wg05**: Forward-Backward scheme described in*J. Opt. B.***7**, 300 (2005). Works for all kinds of target operators, can be used with all kinds of filters, and allows a fixed fluence. The price is a rather unstable convergence. If the restrictions set by the filter and fluence are reasonable, a good overlap can be expected within 20 iterations. No monotonic convergence.**oct_mt03**: Basically an improved and generalized scheme. Comparable to ZBR98/ZR98. See [Y. Maday and G. Turinici,*J. Chem. Phys.***118**, 8191 (2003)].**oct_krotov**: The procedure reported in [D. Tannor, V. Kazakov and V. Orlov, in*Time-Dependent Quantum Molecular Dynamics*, edited by J. Broeckhove and L. Lathouweres (Plenum, New York, 1992), pp. 347-360].**oct_straight_iteration**: Straight iteration: one forward and one backward propagation is performed at each iteration, both with the same control field. An output field is calculated with the resulting wavefunctions.**oct_cg**: Conjugate-gradients, as implemented in the GNU GSL library. In particular, the Fletcher-Reeves version.**oct_bfgs**: The methods use the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. Also, it calls the GNU GSL library version of the algorithm. It is a quasi-Newton method which builds up an approximation to the second derivatives of the function using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region. We have chosen to implement the "bfgs2" version, as GSL calls it, which is supposed to be the most efficient version available, and a faithful implementation of the line minimization scheme described in "Practical Methods of Optimization", (Fletcher), Algorithms 2.6.2 and 2.6.4.**oct_direct**: This is a "direct" optimization scheme. This means that we do not make use of the "usual" QOCT equations (backward-forward propagations, etc), but we use some gradient-free maximization algorithm for the function that we want to optimize. In this case, the maximization algorithm is the Nelder-Mead algorithm as implemeted in the GSL. The function values are obtained by successive forward propagations.

**OCTSpatialCurrWeight**
*Section*: Calculation Modes::Optimal Control
*Type*: block

Can be seen as a position-dependent `OCTCurrentWeight`. Consequently, it
weights contribution of current \(j\) to its functional \(J1_c[j]\) according to the position in space.
For example, `oct_curr_square` thus becomes
\(J1_c[j] = {\tt OCTCurrentWeight} \int{\left| j(r) \right|^2 {\tt OCTSpatialCurrWeight}(r) dr}\).

It is defined as `OCTSpatialCurrWeight`\((r) = g(x) g(y) g(z)\), where
\(g(x) = \sum_{i} 1/(1+e^{-{\tt fact} (x-{\tt startpoint}_i)}) - 1/(1+e^{-{\tt fact} (x-{\tt endpoint}_i)})\).
If not specified, \(g(x) = 1\).

Each \(g(x)\) is represented by one line of the block that has the following form

`%OCTSpatialCurrWeight
dimension | fact | startpoint_1 | endpoint_1 | startpoint_2 | endpoint_2 |...
%`

There are no restrictions on the number of lines, nor on the number of pairs of start- and endpoints. Attention:

**OCTStartIterCurrTg**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: 0

Allows for a time-dependent target for the current without defining it for the total
time-interval of the simulation.
Thus it can be switched on at the iteration desired, `OCTStartIterCurrTg` >= 0
and `OCTStartIterCurrTg` < `TDMaximumIter`.
Tip: If you would like to specify a real time for switching
the functional on rather than the number of steps, just use something
like:
`OCTStartIterCurrTg` = 100.0 / `TDTimeStep`.

**OCTTargetDensity**
*Section*: Calculation Modes::Optimal Control
*Type*: string

If `OCTTargetOperator = oct_tg_density`, then one must supply the target density
that should be searched for. This one can do by supplying a string through
the variable `OCTTargetDensity`. Alternately, give the special string `"OCTTargetDensityFromState"`
to specify the expression via the block `OCTTargetDensityFromState`.

**OCTTargetDensityFromState**
*Section*: Calculation Modes::Optimal Control
*Type*: block
*Default*: no

If `OCTTargetOperator = oct_tg_density`, and `OCTTargetDensity = "OCTTargetDensityFromState"`,
you must specify a `OCTTargetDensityState` block, in order to specify which linear
combination of the states present in `restart/gs` is used to
create the target density.

The syntax is the same as the `TransformStates` block.

**OCTTargetOperator**
*Section*: Calculation Modes::Optimal Control
*Type*: integer
*Default*: oct_tg_gstransformation

The variable `OCTTargetOperator` prescribes which kind of target functional is
to be used.
*Options*:

**oct_tg_velocity**: (Experimental) The target is a function of the velocities of the nuclei at the end of the influence of the external field, defined by`OCTVelocityTarget`**oct_tg_hhgnew**: (Experimental) The target is the optimization of the HHG yield. You must supply the`OCTHarmonicWeight`string. It attempts to optimize the integral of the harmonic spectrum multiplied by some user-defined weight function.**oct_tg_classical**: (Experimental)**oct_tg_spin**: (Experimental)**oct_tg_groundstate**: The target operator is a projection operator on the ground state,*i.e.*the objective is to populate the ground state as much as possible.**oct_tg_excited**: (Experimental) The target operator is an "excited state". This means that the target operator is a linear combination of Slater determinants, each one formed by replacing in the ground-state Slater determinant one occupied state with one excited state (*i.e.*"single excitations"). The description of which excitations are used, and with which weights, should be given in a file called`oct-excited-state-target`. See the documentation of subroutine`excited_states_init`in the source code in order to use this feature.**oct_tg_gstransformation**: The target operator is a projection operator on a transformation of the ground-state orbitals defined by the block`OCTTargetTransformStates`.**oct_tg_userdefined**: (Experimental) Allows to define target state by using`OCTTargetUserdefined`.**oct_tg_jdensity**: (Experimental)**oct_tg_local**: (Experimental) The target operator is a local operator.**oct_tg_td_local**: (Experimental) The target operator is a time-dependent local operator.**oct_tg_exclude_state**: (Experimental) Target operator is the projection onto the complement of a given state, given by the block`OCTTargetTransformStates`. This means that the target operator is the unity operator minus the projector onto that state.**oct_tg_hhg**: (Experimental) The target is the optimization of the HHG yield. You must supply the`OCTOptimizeHarmonicSpectrum`block, and it attempts to optimize the maximum of the spectrum around each harmonic peak. You may use only one of the gradient-less optimization schemes.

**OCTTargetSpin**
*Section*: Calculation Modes::Optimal Control
*Type*: block

(Experimental) Specify the targeted spin as a 3-component vector. It will be normalized.

**OCTTargetTransformStates**
*Section*: Calculation Modes::Optimal Control
*Type*: block
*Default*: no

If `OCTTargetOperator = oct_tg_gstransformation`, you must specify a
`OCTTargetTransformStates` block, in order to specify which linear
combination of the states present in `restart/gs` is used to
create the target state.

The syntax is the same as the `TransformStates` block.

**OCTTargetUserdefined**
*Section*: Calculation Modes::Optimal Control
*Type*: block

Define a target state. Syntax follows the one of the `UserDefinedStates` block.
Example:

`%OCTTargetUserdefined
1 | 1 | 1 | "exp(-r^2)*exp(-i*0.2*x)"
%`

**OCTTdTarget**
*Section*: Calculation Modes::Optimal Control
*Type*: block

(Experimental) If `OCTTargetOperator = oct_tg_td_local`, then you must supply
a OCTTdTarget block. The block should only contain one element, a string cotaining the
definition of the time-dependent local target, *i.e.* a function of x,y,z and t that
is to be maximized along the evolution.

**OCTVelocityDerivatives**
*Section*: Calculation Modes::Optimal Control
*Type*: block

If `OCTTargetOperator = oct_tg_velocity`, and
`OCTScheme = oct_cg` or `OCTScheme = oct_bfgs`
then you must supply the target in terms of the ionic velocities AND
the derivatives of the target with respect to the ionic velocity components.
The derivatives are supplied via strings through the block
`OCTVelocityDerivatives`.
Each velocity component is supplied by `"v[n_atom,vec_comp]"`,
while `n_atom` is the atom number, corresponding to the
`Coordinates` block, and `vec_comp` is the corresponding
vector component of the velocity. The first line of the
`OCTVelocityDerivatives` block contains the derivatives
with respect to `v[1,*]`, the second with respect to `v[2,*]` and so
on. The first column contains all derivatives with respect `v[*,1]`,
the second with respect to `v[*,2]` and the third w.r.t. `v[*,3]`.
As an example, we show the `OCTVelocityDerivatives` block
corresponding to the target shown in the `OCTVelocityTarget`
help section:

`%OCTVelocityDerivatives
" 2*(v[1,1]-v[2,1])" | " 2*(v[1,2]-v[2,2])" | " 2*(v[1,3]-v[2,3])"
"-2*(v[1,1]-v[2,1])" | "-2*(v[1,2]-v[2,2])" | "-2*(v[1,3]-v[2,3])"
%`

**OCTVelocityTarget**
*Section*: Calculation Modes::Optimal Control
*Type*: block

If `OCTTargetOperator = oct_tg_velocity`, then one must supply the
target to optimize in terms of the ionic velocities. This is done by
supplying a string through the block `OCTVelocityTarget`.
Each velocity component is supplied by `"v[n_atom,vec_comp]"`,
where `n_atom` is the atom number, corresponding to the
`Coordinates` block, and `vec_comp` is the corresponding
vector component of the velocity. The target string can be
supplied by using several lines in this block.
As an example, the following target can be used to maximize the
velocity difference between atom 1 and 2 (in a 3D system):

`%OCTVelocityTarget
"(v[1,1]-v[2,1])^2 + (v[1,2]-v[2,2])^2 + "
"(v[1,3]-v[2,3])^2"
%`

**MaximumIter**
*Section*: Calculation Modes::Unoccupied States
*Type*: integer
*Default*: 50

Maximum number of eigensolver iterations. The code will stop even if convergence
has not been achieved. -1 means unlimited. 0 means just do LCAO or read from
restart, and stop.

**UnoccShowOccStates**
*Section*: Calculation Modes::Unoccupied States
*Type*: logical
*Default*: false

If true, the convergence for the occupied states will be shown too in the output.
This is useful for testing, or if the occupied states fail to converge.
It will be enabled automatically if only occupied states are being calculated.