O

Name Occupations
Section States
Type block

The occupation numbers of the orbitals can be fixed through the use of this variable. For example:

%Occupations
  2 | 2 | 2 | 2 | 2
%

would fix the occupations of the five states to 2. There can be at most as many columns as states in the calculation. If there are fewer columns than states, then the code will assume that the user is indicating the occupations of the uppermost states where all lower states have full occupation (i.e. 2 for spin-unpolarized calculations, 1 otherwise) and all higher states have zero occupation. The first column will be taken to refer to the lowest state such that the occupations would be consistent with the correct total charge. For example, if there are 8 electrons and 10 states (from ExtraStates = 6), then an abbreviated specification

%Occupations
  1 | 0 | 1
%

would be equivalent to a full specification

%Occupations
  2 | 2 | 2 | 1 | 0 | 1 | 0 | 0 | 0 | 0
%

This is an example of use for constrained density-functional theory, crudely emulating a HOMO->LUMO+1 optical excitation. The number of rows should be equal to the number of k-points times the number of spins. For example, for a finite system with SpinComponents == spin_polarized, this block should contain two lines, one for each spin channel. All rows must have the same number of columns.

The Occupations block is useful for the ground state of highly symmetric small systems (like an open-shell atom), to fix the occupation numbers of degenerate states in order to help octopus to converge. This is to be used in conjuction with ExtraStates. For example, to calculate the carbon atom, one would do:

ExtraStates = 2
%Occupations
  2 | 2/3 | 2/3 | 2/3
%

If you want the calculation to be spin-polarized (which makes more sense), you could do:

ExtraStates = 2
%Occupations
   2/3 | 2/3 | 2/3
   0 | 0 | 0
%

Note that in this case the first state is absent, the code will calculate four states (two because there are four electrons, plus two because ExtraStates = 2), and since it finds only three columns, it will occupy the first state with one electron for each of the spin options.

If the sum of occupations is not equal to the total charge set by ExcessCharge, an error message is printed. If FromScratch = no and RestartFixedOccupations = yes, this block will be ignored.

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

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

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

Name 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)
  

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

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

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

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

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

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

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

Name 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).

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

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

Name 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) frequency_filter: Specifies a spectral filter.

(ii) time_filter: DISABLED IN THIS VERSION.

Example:

%OCTFilter
  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.

Options:

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

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

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

Name 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)$.

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

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

Name 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(-i0.2*x)"
%

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

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

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

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

Name 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].

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

Name 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
%

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

Name 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].

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

Name OCTScheme
Section Calculation Modes::Optimal Control
Type integer
Default oct_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_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. The seed for the random number generator can be modified by setting GSL_RNG_SEED environment variable.
  
• 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 implemented in the GSL. The function values are obtained by successive forward propagations. The seed for the random number generator can be modified by setting GSL_RNG_SEED environment variable.
  
• 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. The seed for the random number generator can be modified by setting GSL_RNG_SEED environment variable.
  
• 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. The seed for the random number generator can be modified by setting GSL_RNG_SEED environment variable.
  

Name 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: startpoint and endpoint have to be supplied pairwise with startpoint < endpoint. dimension > 0 is integer, fact is float.

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

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

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

Name 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_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_elec_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.
  
• 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)
  

Name OCTTargetSpin
Section Calculation Modes::Optimal Control
Type block

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

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

Name 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(-i0.2*x)"
%

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

Name 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])"
%

Name 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"
%

Name OEPLevel
Section Hamiltonian::XC
Type integer
Default oep_kli

At what level shall Octopus handle the optimized effective potential (OEP) equation.

Options:

• oep_none:
  Do not solve OEP equation.
  
• oep_kli:
  Krieger-Li-Iafrate (KLI) approximation. Ref: JB Krieger, Y Li, GJ Iafrate, Phys. Lett. A 146, 256 (1990).
  
• oep_full:
  (Experimental) Full solution of OEP equation using the Sternheimer approach. The linear solver will be controlled by the variables in section Linear Response::Solver, and the iterations for OEP by Linear Response::SCF in LR calculations and variable OEPMixing. Note that default for LRMaximumIter is set to 10. Ref: S. Kuemmel and J. Perdew, Phys. Rev. Lett. 90, 043004 (2003).
  

Name OEPMixing
Section Hamiltonian::XC
Type float
Default 1.0

The linear mixing factor used to solve the Sternheimer equation in the full OEP procedure.

Name OEPMixingScheme
Section Hamiltonian::XC
Type integer
Default oep_mixing_scheme_const

Different Mixing Schemes are possible

Options:

• oep_mixing_scheme_const:
  Use a constant Reference: S. Kuemmel and J. Perdew, Phys. Rev. Lett. 90, 4, 043004 (2003)
  
• oep_mixing_scheme_bb:
  Use the Barzilai-Borwein (BB) Method Reference: T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos, Phys. Rev. B 85, 235126 (2012)
  
• oep_mixing_scheme_dens:
  Use the inverse of the electron density Reference: S. Kuemmel and J. Perdew, Phys. Rev. B 68, 035103 (2003)
  

Name OEPRemoveElectron
Section Hamiltonian::XC
Type logical
Default .false.

Remove electron-electron interaction in OEP-Photon calculations

Name OnlyUserDefinedInitialStates
Section States
Type logical
Default no

If true, then only user-defined states from the block UserDefinedStates will be used as initial states for a time-propagation. No attempt is made to load ground-state orbitals from a previous ground-state run.

Name OperateAccel
Section Execution::Optimization
Type integer
Default map

This variable selects the subroutine used to apply non-local operators over the grid when an accelerator device is used.

Options:

• invmap:
  The standard implementation ported to OpenCL.
  
• map:
  A different version, more suitable for GPUs.
  

Name OperateComplex
Section Execution::Optimization
Type integer
Default optimized

This variable selects the subroutine used to apply non-local operators over the grid for complex functions.

Options:

• fortran:
  The standard Fortran function.
  
• optimized:
  This version is optimized using vector primitives (if available).
  

Name OperateDouble
Section Execution::Optimization
Type integer
Default optimized

This variable selects the subroutine used to apply non-local operators over the grid for real functions.

Options:

• fortran:
  The standard Fortran function.
  
• optimized:
  This version is optimized using vector primitives (if available).
  

Name OptimizeChebyshevFilterDegree
Section SCF::Eigensolver
Type logical
Default yes

Used by the Chebyshev filter only. Octopus generates a best-estimate for the degree of the Chebyshev polynomial used to filter the subspace. A separate estimate is generated for each state block, per SCF iteration. Note that if running parallelism over states, the block/batch containing the largest eigenstates will converge more slowly and will therefore use a larger degree relative to all other batches. One should therefore avoid setting "ChebyshevFilterDegree" to an excessive value (for example > 50). For more details regarding how the degree is estimated, one can refer to Section 5.4 in "Computer Physics Communications 187 (2015) 98–105" (http://dx.doi.org/10.1016/j.cpc.2014.10.015).

Name Output
Section Output
Type block
Default none

Specifies what to print. Each output must be in a separate row. Optionally individual output formats and output intervals can be defined for each row or they can be read separately from OutputFormat and OutputInterval variables in the input file. The output files are written at the end of the run into the output directory for the relevant kind of run (e.g. static for CalculationMode = gs). Time-dependent simulations print only per iteration, including always the last. The frequency of output per iteration (available for CalculationMode = gs, unocc, td, and opt_control) is set by OutputInterval and the directory is set by OutputIterDir. For linear-response run modes, the derivatives of many quantities can be printed, as listed in the options below. Indices in the filename are labelled as follows: sp = spin (or spinor component), k = k-point, st = state/band. There is no tag for directions, given as a letter. The perturbation direction is always the last direction for linear-response quantities, and a following +/- indicates the sign of the frequency.

Example (minimal):

%Output
  density
  potential
%

Example (with OutputFormat):

%Output
  density | cube + axis_z
  potential | cube
%

Example (with OutputFormat, incomplete):

%Output
  density | cube + axis_z
  potential
%

Example (tagged):

%Output
  density | "output_format" | cube + axis_z | "output_interval" | 50
  potential | "output_format" | cube | "output_interval" | 20
%

Example (tagged, incomplete):

%Output
  density | "output_format" | cube + axis_z
  potential | "output_interval" | 20
%
Missing information for the incomplete blocks will be parsed form the out-of-block definitions. It is also possible to mix the order of columns in the tagged format. See OutputFormat, and OutputInterval.

Options:

• potential:
  Outputs Kohn-Sham potential, separated by parts. File names are v0 for the local part of the ionic potential, vc for the classical potential (if it exists), vh for the Hartree potential, vks for the local part of the Kohn-Sham potential, and vxc- for the exchange-correlation potentials. For vks and vxc, a suffix for spin is added in the spin-polarized case.
  
• density:
  Outputs density. The output file is called density-, or lr_density- in linear response.
  
• wfs:
  Outputs wavefunctions. Which wavefunctions are to be printed is specified by the variable OutputWfsNumber -- see below. The output file is called wf-, or lr_wf- in linear response.
  
• wfs_sqmod:
  Outputs modulus squared of the wavefunctions. The output file is called sqm-wf-. For linear response, the filename is sqm_lr_wf-.
  
• geometry:
  Outputs file containing the coordinates of the atoms treated within quantum mechanics. If OutputFormat = xyz, the file is called geometry.xyz; a file crystal.xyz is written with a supercell geometry if the system is periodic; if point charges were defined in the PDB file (see PDBCoordinates), they will be output in the file geometry_classical.xyz. If OutputFormat = xcrysden, a file called geometry.xsf is written.
  
• current:
  Outputs the total current density. The output file is called current-. For linear response, the filename is lr_current-.
  
• ELF:
  Outputs electron localization function (ELF). The output file is called elf-, or lr_elf- in linear response, in which case the associated function D is also written, as lr_elf_D-. Only in 2D and 3D.
  
• ELF_basins:
  Outputs basins of attraction of the ELF. The output file is called elf_rs_basins.info. Only in 2D and 3D.
  
• Bader:
  Outputs Laplacian of the density which shows lone pairs, bonded charge concentrations and regions subject to electrophilic or nucleophilic attack. See RF Bader, Atoms in Molecules: A Quantum Theory (Oxford Univ. Press, Oxford, 1990).
  
• el_pressure:
  Outputs electronic pressure. See Tao, Vignale, and Tokatly, Phys Rev Lett 100, 206405 (2008).
  
• matrix_elements:
  Outputs a series of matrix elements of the Kohn-Sham states. What is output can be controlled by the OutputMatrixElements variable.
  
• pol_density:
  Outputs dipole-moment density dipole_density-, or polarizability density alpha_density- in linear response. If ResponseMethod = finite_differences, the hyperpolarizability density beta_density- is also printed.
  
• mesh_r:
  Outputs values of the coordinates over the grid. Files will be called mesh_r- followed by the direction.
  
• kinetic_energy_density:
  Outputs kinetic-energy density, defined as:

  $\tau_\sigma(\vec{r}) = \sum_{i=1}^{N_\sigma} \left| \vec{\nabla} \phi_{i\sigma}(\vec{r}) \right|^2,. $

  The index $\sigma$ is the spin index for the spin-polarized case, or if you are using spinors. For spin-unpolarized calculations, you get the total kinetic-energy density. The previous expression assumes full or null occupations. If fractional occupation numbers, each term in the sum is weighted by the occupation. Also, if we are working with an infinite system, all k-points are summed up, with their corresponding weights. The files will be called tau-sp1 and tau-sp2, if the spin-resolved kinetic energy density is produced (runs in spin-polarized and spinors mode), or only tau if the run is in spin-unpolarized mode.
  

• dos:
  Outputs density of states. See DOSEnergyMax, DOSEnergyMin, DOSEnergyPoints, and DOSGamma.
  
• tpa:
  Outputs transition-potential approximation (TPA) matrix elements, using $\vec{q}$-vector specified by MomentumTransfer.
  
• forces:
  Outputs file forces.xsf containing structure and forces on the atoms as a vector associated with each atom, which can be visualized with XCrySDen.
  
• wfs_fourier:
  (Experimental) Outputs wavefunctions in Fourier space. This is only implemented for the ETSF file format output. The file will be called wfs-pw-etsf.nc.
  
• xc_density:
  Outputs the XC density, which is the charge density that generates the XC potential. (This is $-1/4\pi$ times the Laplacian of the XC potential). The files are called nxc.
  
• PES_wfs:
  Outputs the photoelectron wavefunctions. The file name is pes_wfs- plus the orbital number.
  
• PES_density:
  Outputs the photolectron density. Output file is pes_dens- plus spin species if spin-polarized calculation is performed.
  
• PES:
  Outputs the time-dependent photoelectron spectrum.
  
• BerkeleyGW:
  Output for a run with BerkeleyGW. See Output::BerkeleyGW for further specification.
  
• delta_perturbation:
  Outputs the "kick", or time-delta perturbation applied to compute optical response in real time.
  
• external_td_potential:
  Outputs the (scalar) time-dependent potential.
  
• mmb_wfs:
  Triggers the ModelMB wavefunctions to be output for each state.
  
• mmb_den:
  Triggers the ModelMB density matrix to be output for each state, and the particles specified by the DensitytoCalc block. Calculates, and outputs, the reduced density matrix. For the moment the trace is made over the second dimension, and the code is limited to 2D. The idea is to model N particles in 1D as an N-dimensional non-interacting problem, then to trace out N-1 coordinates.
  
• potential_gradient:
  Prints the gradient of the potential.
  
• energy_density:
  Outputs the total energy density to a file called energy_density.
  
• heat_current:
  Outputs the total heat current density. The output file is called heat_current-.
  
• photon_correlator:
  Outputs the electron-photon correlation function. The output file is called photon_correlator.
  
• J_flow:
  todo: document J_flow option!
  
• current_kpt:
  Outputs the current density resolved in momentum space. The output file is called current_kpt-.
  
• density_kpt:
  Outputs the electronic density resolved in momentum space.
  
• occ_matrices:
  Only for DFT+U calculations. Outputs the occupation matrices of DFT+U
  
• effectiveU:
  Only for DFT+U calculations. Outputs the value of the effectiveU for each atoms
  
• magnetization:
  Only for DFT+U calculations. Outputs file containing structure and magnetization of the localized subspace on the atoms as a vector associated with each atom, which can be visualized. For the moment, it only works if a +U is added on one type of orbital per atom.
  
• local_orbitals:
  Only for DFT+U calculations. Outputs the localized orbitals that form the correlated subspace
  
• kanamoriU:
  Only for DFT+U calculations. Outputs the Kanamori interaction parameters U, U`, and J. These parameters are not determined self-consistently, but are taken from the occupation matrices and Coulomb integrals comming from a standard +U calculation.
  
• xc_torque:
  Outputs the exchange-correlation torque. Only for the spinor case and in the 3D case.
  
• eigenval_kpt:
  Outputs the eigenvalues resolved in momentum space, with one file for each band.
  
• stress:
  Outputs the stress tensor and each of its contributing terms
  
• current_dia:
  Outputs the diamagnetic current density from a non-uniform vector potential. The output file is called current_dia-.
  
• jdos:
  Outputs the joint density of states. The same variables as for the regular DOS are used to control the energies and broadening.
  

Name OutputDuringSCF
Section Output
Type logical
Default no

During gs and unocc runs, if this variable is set to yes, output will be written after every OutputInterval iterations.

Name OutputFormat
Section Output
Type flag
Default 0

Describes the format of the output files. This variable can also be defined inside the Output block. See Output. Example: axis_x + plane_x + dx

Options:

• axis_x:
  The values of the function on the x axis are printed. The string .y=0,z=0 is appended to previous file names.
  
• axis_y:
  The values of the function on the y axis are printed. The string .x=0,z=0 is appended to previous file names.
  
• axis_z:
  The values of the function on the z axis are printed. The string .x=0,y=0 is appended to previous file names.
  
• plane_x:
  A plane slice at x = 0 is printed. The string .x=0 is appended to previous file names.
  
• plane_y:
  A plane slice at y = 0 is printed. The string .y=0 is appended to previous file names.
  
• plane_z:
  A plane slice at z = 0 is printed. The string .z=0 is appended to previous file names.
  
• dx:
  For printing three-dimensional information, the open-source program visualization tool OpenDX can be used. The string .dx is appended to previous file names. Available only in 3D.
  
• netcdf:
  Outputs in NetCDF format. This file can then be read, for example, by OpenDX. The string .ncdf is appended to previous file names. Requires the NetCDF library. Only writes the real part of complex functions.
  
• mesh_index:
  Generates output files of a given quantity (density, wavefunctions, ...) which include the internal numbering of mesh points. Since this mode produces large datafiles this is only useful for small meshes and debugging purposes. The output can also be used to display the mesh directly. A Gnuplot script for mesh visualization can be found under PREFIX/share/octopus/util/display_mesh_index.gp.
  
• xcrysden:
  A format for printing structures and three-dimensional information, which can be visualized by the free open-source program XCrySDen and others. The string .xsf is appended to previous file names. Note that lattice vectors and coordinates are as specified by UnitsOutput. Available in 2D and 3D.
  
• matlab:
  In combination with plane_x, plane_y and plane_z, this option produces output files which are suitable for 2D Matlab functions like mesh(), surf(), or waterfall(). To load these files into Matlab you can use, e.g. >> density = load('static/density-1.x=0.matlab.abs'); >> mesh(density);
  
• meshgrid:
  Outputs in Matlab mode the internal mesh in a format similar to >> [x,y] = meshgrid(-2:.2:2,-1:.15:1) The x meshgrid is contained in a file *.meshgrid.x and the y-grid can be found in *.meshgrid.y.
  
• boundary_points:
  This option includes the output of the mesh enlargement. Default is without. Supported only by binary, axis, plane, mesh_index, and matlab formats. Not all types of Output will have this information available. Not supported when parallel in domains.
  
• binary:
  Plain binary, new format.
  
• etsf:
  ETSF file format. Requires the ETSF_IO library. Applies only to Output = density, geometry, wfs, and/or wfs_fourier.
  
• xyz:
  Geometry will be output in XYZ format. Does not affect other outputs.
  
• cube:
  Generates output in the cube file format. Available only in 3D. Only writes the real part of complex functions. This output format always uses atomic units.
  
• bild:
  Generates output in BILD format.
  
• vtk:
  Generates output in VTK legacy format.
  
• integrate_xy:
  Integrates the function in the x-y plane and the result on the z axis is printed.
  
• integrate_xz:
  Integrates the function in the x-z plane and the result on the y axis is printed
  
• integrate_yz:
  Integrates the function in the y-z plane and the result on the x axis is printed
  
• ascii:
  Plain text format regardless of dimensionality. For the moment only employed by the oct-phototoelectron_spectrum post-processing utility.
  

Name OutputInterval
Section Output
Type integer
Default 50

The output requested by variable Output is written to the directory OutputIterDir when the iteration number is a multiple of the OutputInterval variable. Subdirectories are named Y.X, where Y is td, scf, or unocc, and X is the iteration number. To use the working directory, specify "." (Output of restart files is instead controlled by RestartWriteInterval.) Must be >= 0. If it is 0, then no output is written. For gs and unocc calculations, OutputDuringSCF must be set too for this output to be produced. This variable can also be defined inside the Output block. See Output.

Name OutputIterDir
Section Output
Type string
Default “output_iter”

The name of the directory where Octopus stores information such as the density, forces, etc. requested by variable Output in the format specified by OutputFormat. This information is written while iterating CalculationMode = gs, unocc, or td, according to OutputInterval, and has nothing to do with the restart information.

Name OutputMatrixElements
Section Output
Type block
Default none

Specifies what matrix elements to print. Enabled only if Output block includes matrix_elements. The output files go into the static directory, except when running a time-dependent simulation, when the directory td.XXXXXXX is used.

Example:

%OutputMatrixElements
  momentum
  ks_multipoles
%

It is possible to specify only compute the matrix elements for some of the states using the variables OutptMEStart and OutputMEEnd.

Options:

• momentum:
  Momentum. Filename: ks_me_momentum.
  
• ang_momentum:
  Dimensionless angular momentum $\vec{r} \times \vec{k}$. Filename: ks_me_angular_momentum.
  
• one_body:
  $\left< i \left| \hat{T} + V_{ext} \right| j \right>$. Not available with states parallelization.
  
• two_body:
  $\left< ij \left| \frac{1}{\left|\vec{r}_1-\vec{r}_2\right|} \right| kl \right>$. Not available with states parallelization. Not available with states parallelization. For periodic system, this is not available for k-point parallelization neither.
  
• two_body_exc_k:
  $\left< n1-k1, n2-k2 \left| \frac{1}{\left|\vec{r}_1-\vec{r}_2\right|} \right| n2-k1 n1-k2 \right>$. Not available with states parallelization. For periodic system, this is not available for k-point parallelization neither.
  
• ks_multipoles:
  See OutputMEMultipoles. Not available with states parallelization.
  
• dipole:
  Prints the dipole matrix elements. Not available with states parallelization. For periodic systems, the intraband terms (dipole matrix elements between degenerated states) are set to zero, and only the absolute value of the dipole matrix element is printed. Not yet supported for spinors.
  

Name OutputMEEnd
Section Output
Type integer
Default 1

Specifies the highest state/band index used to compute the matrix element. So far, this is only used for dipole matrix elements.

Name OutputMEMultipoles
Section Output
Type integer
Default 1

This variable decides which multipole moments are printed out for OutputMatrixElements = ks_multipoles:

In 3D, if, for example, OutputMEMultipoles = 1, then the program will print three files, ks_me_multipoles.x (x=1,2,3), containing respectively the (1,-1), (1,0) and (1,1) multipole matrix elements between Kohn-Sham states.

In 2D, this variable is ignored: it will always print two files, ks_me_multipoles.i (i=1,2), containing the $x$ and $y$ dipole matrix elements.

In 1D, if, for example, OutputMEMultipoles = 2, the program will print two files, containing the $x$ and $x^2$ matrix elements between Kohn-Sham states.

Name OutputMEStart
Section Output
Type integer
Default 1

Specifies the state/band index for starting to compute the matrix element. So far, this is only used for dipole matrix elements.

Name OutputWfsNumber
Section Output
Type string
Default all states

Which wavefunctions to print, in list form: i.e., "1-5" to print the first five states, "2,3" to print the second and the third state, etc. If more states are specified than available, extra ones will be ignored.