The fitting of the parameters is done using the VCHFIT program. This reads an input file, specified as the first argument
vchfitXX input
or
vchfitXX name.inp
where XX is the version number, e.g. vchfit83 for package version 8.3, and where name (or name.inp) denotes the input file. The input format uses, like all the MCTDH package input files, keywords that are for the most part free format and case insensitive. See MCTDH input file structure

for further information on the general use of keywords, noting that there are no sections in the VCHAM input files. The input file ends with the keyword end-input

Input file (*.inp)

The input file contains a set of keywords to control exactly how the fit is made. The following list is of the groups of keywords that are used

Defining the Data Set

The basic ingredient for a fit is a set of data - energies and derivatives at various points. These are contained in a *.info file set up previously (see setting up a database )

Defining the Data set
Keyword Argument Description
infofile S S is the name of the .info file

Defining the Model

The vibronic coupling model Hamiltonian contains sets of parameters which can be defined as on-diagonal and off-diagonal. These sets can be further divided according to the power of the coordinate: linear, second-order, etc. In the following documentation the sets of parameters are described using the nomenclature in the following table:
Sets of Parameters
Name Description
E (s,s1) Constant for states s,s1
kappa (m,s) On-diagonal linear parameters for mode m state s
lambda (m,s,s1) Off-diagonal linear parameters for mode m connecting states s and s1
gamma (m,m1,s) On-diagonal second-order parameters for modes m,m1 and state s
mu (m,m1,s,s1) Off-diagonal second-order parameters for modes m,m1 connecting states s and s1
general (k,m,s) The expansion coefficient k for a general potential for mode m and state s
cpgeneral (k,m,s,s1) The expansion coefficient k for a general diabatic coupling function for mode m and states s and s1

Which sets of parameters are included in the fit can be controlled by the following keywords:

Defining the Model
Keyword Description
linear_mod linear model - Constants + kappa + lambda
adiab2_mod Adiabatic 2nd order model - Constants + kappa + lambda + gamma
diab2_mod Diabatic 2nd order model - Constants + kappa + lambda + gamma + mu

Further control is possible by selecting sets of the parameters for fitting:
Choosing the parameters to fit
Keyword Description
fullfit All parameters are included
linfit All constants and linear terms included
constfit Diagonal constants included
constall All constants included
kappafit On-diagonal linear terms
lambdafit Off-diagonal linear terms
gammafit On-digonal second-order terms (quadratic and bi-linear)
gammafit_diag On-diagonal quadratic terms
mufit Off-diagonal second-order terms
iotafit Only on-diagonal quadratic-linear terms
iota2fit Only off-diagonal quadratic-linear terms
highfit = S On-diagonal higher order terms. S = 3 in addition to iotafit completes 3rd order. S=4 adds fourth-order
diabfit Coefficients for diabatic functions
read_mask The symmetry "mask" defining which parameters are to be optimised will be read from a .vcham file specified by the guess = S keyword
Thus if one chooses diab2_mod linfit all parameters up to second-order both on- and off-diagonal are included in the Hamiltonian, but only parameters up to linear will be optimised. In this way all parameters are written to the output files and can be used in later stages.

More detailed control is possible using keywords described below, or by editing and using a .vcham file to set up the initial guess.

It is also possible to fit only a set of modes.

Selecting Modes
Keyword Description
fit_all Fit all modes (default)
fit_select = m [,m1 [, m2 ....]] Fit only parameters for modes m [and m1 etc.]

Defining the Non-zero Parameters

If a system has symmetry, many parameters in the model Hamiltonian are zero from group theory arguments. The program reads information defining the system from the *.info file, e.g. the number of electronic states and vibrational modes in the system. It also reads the symmetry labels for the modes. There are two options for the determination of parameters that are necessarily zero by symmetry:

(i) Specification of the symmetry labels of the modes that couple given pairs of states

The following keywords use the symmetry labels to define which parameters are non-zero by symmetry. In each case label is the symmetry label as written in the *.info file.

Defining the Non-zero Parameters
Keyword Description
symmode = label "label" is the symmetry label for the totally symmetric modes
cpmode = s s1 label "label" is the symmetry label for the modes coupling states s and s1
cpmode2 = s s1 label label1 "label" and "label1" are the symmetries for the pairs of modes coupling states s and s1
jtmode = s label "label" is the symmetry label for any Jahn-Teller active modes for state s
jtstates = s, s1 [,s2, ...] States s and s1 are a degenerate pair with Jahn-Teller active modes defined
using jtmode and coupled by the mode defined by cpmode(s,s1). Higher
order JT than double degenrate is possible if s2 etc. defined.

The symmode label defines the modes that have non-zero on-diagonal linear parameters (all modes have non-zero on-diagonal second-order parameters). The cpmode(s,s1) label is the product of the symmetries of states s and s1. Modes with this symmetry have non-zero linear off-diagonal (coupling) parameters. Two modes whose product has this symmetry will have non-zero second-order off-diagonal terms (not implemented fully). If the system has a Jahn-Teller intersection, then in addition to the totally symmetric modes, modes with symmetry label jtmode(s) can appear as linear parameters on the diagonal of state s.

(ii) Specification of state symmetries

By calculating the number of times the totally symmetric irreducible representation enters into the decomposition of the direct product representations generated by the integrands of the parameters of the vibronic coupling Hamiltonian, the VCHFIT program can determine automatically whether or not a parameter is zero by symmetry. To do so, the point group of the molecule and the symmetry species of the states have to be specified in the input. All Albelian point groups are supported, i.e., the groups Cs, Ci, C2, D2, C2v, C2h, and D2h.

Only Abelian point group symmetry may be used, so for molecules with degenerate symmetry, an Abelian subgroup must be used.

The point group of the molecule is given using the syntax

Specifying the point group
Keyword Description
point_group = PG "PG" is the point group of the molecule

The symmetry species of the states are given between the lines

state_symmetry
and
end-state_symmetry

using the following syntax:

Specifying state symmetries
Keyword Description
s label "s" is the state number, "label" is the symmetry label

By way of example, for a molecule belonging to the C2v point group and three states, the following could possibly be specified:

point_group = C2v

state_symmetry
1 A1
2 B1
3 B2
end-state_symmetry

If any other parameters need to be added to the calculation they are set in a list between the lines
add_parameters-section
and
end-add_parameters-section
(addpar can be used in place of the rather long add_parameters). This can be used, e.g. to add a mu parameter in a system for which the bilinear symmetry rules are not known (in this case diab_mod2 must be selected), or to explicitely add a term of any power, even higher than 3rd.

The parameters are added using the following syntax:

Explicitely adding parameters
Keyword Description
E s s1 = S add a constant
kappa m s = S add an on-diagonal linear parameter
lambda m s s1 = S add an off-diagonal linear parameter
gamma m m1 s = S add an on-diagonal second-order parameter
mu m m1 s s1 = S add an off-diagonal second-order parameter
explicit s s1 m^p m1^p1 ... add explicit term where m, m1, .... define coordinates raised to the powers p, p1,.... (see examples below).
S = optimize Added parameter will be optimised
S = freeze Added parameter will be added but kept frozen
Examples of explicit terms that can be added are addpar-section
explicit 1 1 1^2 3^4
explicit 2 2 1^2 2 3^2
explicit 2 3 8^3
explicit 2 3 8^4
end-addpar-section
The first line adds a quadratic-quartic term coupling modes 1 and 3, on-diagonal in state 1.
The second line adds a quadratic-linear-quadratic terms between modes 1, 2 and 3 on-diagonal in state 2. Note that the linear power can be omitted.
The third and fourth lines add a cubic and quartic term to the inter-state coupling between states 2 and 3 along mode 8.

Defining the Zero-Order Diabatic Potential

By default, the zero-order diabatic potential is a harmonic oscillator function for each mode with the frequency taken from the .info file. A different function can be selected from the following table and written between the keywords

diabatic_function
and
end-define

Zero-Order Diabatic Potential Functions
Keyword Description
m s quartic Quartic potential
m s morse Morse oscillator
m s antimorse Anti-Morse oscillator
m s exp Exponential (unbound) potential
m s avoided_crossing An avoided crossing model potential between an upper exponential and a lower Morse potential, with constant coupling between the two
m s general morse A Morse oscillator that is defined without reference to the frequency
m s morse_lorentzian Sum of a Morse oscillator and a Lorentzian
m s morse_gaussian Sum of a Morse oscillator and a Gaussian
m s avoided_crossing2 An avoided crossing model potential between an upper exponential and a lower Morse potential, with coupling between the two that is a function of the mode m and suitable for use with a dissociative coordinate
m s avoided_crossing_quartic An avoided crossing model potential between an upper quartic potential and a lower quartic potential, with linear coupling between the two

The definitions of the potential functions are here .

Defining the Diabatic Coupling Function

By default, the off-diagonal elements of the diabatic potential are expanded as a Taylor series. It is, however, possible to replace the Taylor expansion along a single mode by a different one-mode function. A different function can be selected from the following table and written between the keywords

coupling_function
and
end-coupling_function

Using the syntax

m s s1 func

where m is the mode number, s and s1 are the state numbers, and func is the name of the function to be used.

Function names and forms are described here.

Switching functions

By default, two modes are correlated by those terms entering into the Taylor expansion of the diabatic potential with respect to both those two modes. Consider the correlation of a bound mode, Qb, and a dissociative mode, Qd. We can anticipate large amplitude motion along Qd, and so the description of the correlation of Qb and Qd via a low-order Taylor expansion should be inadequate. In such cases, it is possible to replace the terms correlating the two modes (bilinear gamma terms, etc.) with appropriate switching functions that have an approproate asymptotic behaviour.

Switching functions are defined between the keywords switching_function
and end-switching_function as follows

switching_function
m parameter func
.
.
.
end-switching_function

Here m is the dissociative mode number, parameter is the name of the parameter associated with the bound mode (specified using the syntax detailed here), and func specifies the switching function to use. For example,

1 kappa 2 3 tanh

would result in the kappa term for mode 2 and state 3 being expanded as a function of mode 1 using the switching function labeled tanh

The switching function labels and their forms are given here .

Selecting Data Sets

The data sets in the .info file can be pruned by using the following keywords :
Selecting Data Sets
Keyword Description
eupper = R Data sets with an energy > R are ignored

Setting Initial Values for the Parameters

The initial values of parameters can be set by reading in the results of a previous fit. They may be also be selectively set to zero.
Setting Initial values for Parameters
Keyword Description
guess = S The initial values are read from the .vcham file S
nonopt_zero All parameters that are not to be optimised are set to zero, irrespective of whether they are zero by symmetry or not.
nozero_zero By default, all parameters zero by symmetry are set to zero before the optimisation is done. This keyword turns off the initial zeroing so that if a value has been read in for one of these parameters it is kept at that value.

The initial value for a parameter may be defined by listing the desired values between the strings

set_parameters-section
and
end-set_parameters-section

VALUE in the following table can be any number and energy units may be used as usual. The description uses the parameter nomenclature from above.

Setting the initial value for a parameter
Keyword Description
E s s1 = VALUE set E(s,s1) to VALUE
kappa m s = VALUE set kappa(m,s) to VALUE
lambda m s s1 = VALUE set lambda(m,s,s1) to VALUE
gamma m m1 s = VALUE set gamma(m,m1,s) to VALUE
mu m m1 s s1 = VALUE set mu(m,m1,s,s1) to VALUE
general k m s = VALUE set general(k,m,s) to VALUE
cpgeneral k m s s1 = VALUE set cpgeneral(k,m,s,s1) to VALUE
switch k m pname = VALUE set the kth coefficient of the switching function defined with respect to mode m for the parameter labeled 'pname' to VALUE, e.g., switch k m kappa m1 s would correspond to the variation of kappa(m1,s) with respect to the mode m.

Constraining the Parameters

Constraints can be set to keep relationships between parameters during the optimisation. To use constraints, the following keyword must be used:
Constraining the parameters
Keyword Description
use_cons use the constraints defined either in input file or vchfit file from which symmetry mask is read

If Jahn-Teller modes have been defined the program automatically tries to set up the appropriate contraints. Other constraints are set in a list between the lines
constraints-section
and
end-constraints-section

The constraints are set using the following syntax:

Constraining the parameters
Keyword Description
E s s1 = CONSTRAINT constrain a constant
kappa m s = CONSTRAINT constrain an on-diagonal linear parameter
lambda m s s1 = CONSTRAINT constrain an off-diagonal linear parameter
gamma m m1 s = CONSTRAINT constrain an on-diagonal second-order parameter
mu m m1 s s1 = CONSTRAINT constrain an off-diagonal second-order parameter
general k m s = CONSTRAINT constrain an expansion coefficient of a general potential

where using the parameter nomenclature above CONSTRAINT can be one of the following:

Possible Constraints
CONSTRAINT Description
freeze keep the parameter fixed
unfreeze allow a previously fixed parameter to be optimised
E s s1 factor constrain to the factor*constant E(s,s1)
kappa m s factor constrain to factor*kappa(m,s)
lambda m s s1 factor constrain to factor*lambda(m,s,s1)
gamma m m1 s factor constrain to factor*gamma(m,m1,s)
mu m m1 s s1 factor constrain to factor*mu(m,m1,s,s1)
general k m s factor constrain to factor*general(k,m,s)
For example to set the on-diagonal linear parameters of mode 1 equal and opposite in two states and equal to the linear off-diagonal parameter of mode 2
Constraints-section
kappa 1 2 = kappa 1 1 -1.0
kappa 1 2 = lambda 2 1 2 1.0
End-Constraints-section

Controlling the Optimisation

Various features of the optimisation can be controlled. In particular, factors can be chosen to weight the points in the database, e.g. by energy (lower energy more important)
Weighting the Points
Keyword Description
equalweight The points are equally weighted (default)
energyweight [ = R ] The points are exponentially weighted according to their energy
wi = exp -R(Ei - E0). If the argument R is not input, R=1.0

The optimisation algorithm may be selected

Controlling the Optimisation
Keyword Description
iter = I I iterations will be made
tol = R The convergence tolerance is set to R
optimiser = S The optimser used will be S where
Keyword Description
conj_grad Conjugate Gradient (default). A simple routine using derivatives.
powell Powell. This uses second derivatives.
simplex Simplex. Requires no derivatives.
genetic Genetic Algorithm. Requires no derivatives, see Genetic algoithm section below
no-ortho The normal mode vectors are not initially orthonormalised but used as read in.
nofit Do not make the fit - only a guess is set up. This differs in effect from setting iter=0 as any constraints set up or modes selected will be ignored.

Genetic Algorithm optimisation

Genetica Algorithm optimisation algorithm uses ideas from evolutionary biology such as inheritance,mutation, selection and crossover in order to reach a soloution. This is very different from the other optimisation techniques which improve on a given guess, and leave the local minima they start in. The Genetic Algorithm creates a simulation of a population representing possible fits, and selects those with the best fit to proceed to the next generation, over . New fits are created through crossover (breeding) between two members of the population and random mutations. For a more detailed description of Genetic Algorithms in general please see this wikipedia entry.

By running this simulation over many generations (typically 100-200) one can be confident that the end fit is in the global rather than a local minimum. As the genetic algorithm is typically time consuming, especially on larger systems, running the algorithm to convergence is not recommended as once a the global minima is reached numerical optimisation (such as those above) is far more efficient.

The following table gives the options available when running Genetic Algorithm Optimisations.

populationThe initial population number, and the number each generation will be reduced to after selection. This number must be even. The maximum number of fits generated is 2 * population * number of mutations. Default is 50, a population below 40 is not recommended and a higher population is only desirable above 250 optimisable parameters.
generations The maximum number of generations the optimisation can run for. Default is 150 but unlike population can vary signifcantly with the number of optimisable parameters.
mutateThe number of population members to be randomly mutated, must be an even number. Default is 6, and a mutation number of around 10% of the population is desirable.
mutation_type = x y This specifies the type of scaling to be used when mutating. As the population approaches convergence the mutation needs to be scaled down otherwise they will be lost from the population within a few generations leading to a lack of diversity in the population. y is the magnitude of the change, taking a value between 1-5 (default 4). As mutations are also scaled by frequency and (as shown below) by generation a direct relationship is complex.
xScaling type
1Stepped, mutation is decreased by a tenth every tenth of the maximum number of generations
2Linear, mutation is decreased linearly starting at y and decreasing to a tenth of y
3None, mutation is not scaled
rngseed = x As computers are incapable of producing true random numbers psuedo random number generators are used. If two runs of the algorithm are made using the same seed numbers the results will be identical, this can be changed by altering the seed x (which may have any value between 1 and 2,147,483,562).

The genetic algorithm stores all potential fits on the hard drive rather than in memory due to the large ammount of data required. These are stored as .vcham files (without the mask portion) with the format xxxx_yyyy.vcham where xxxx is the population number and yyyy is the generation number. Only the last ten generations are stored to conserve hard drive space, this can be changed by searching for ' deleteing all but the last ten generations' in the source code.

Output files

Various output files are produced. These include

Output files
file name Description
name.vcham contains the data from the fit
name.log the log file
name.par the parameters in a Quantics .par file format. To write this file the keyword quantics_par must be set
name.op the operator in a Quantics .op file format. to write this file the keyord quantics_op must be set
name.inp a simple Quantics .inp file. to write this file the keyord quantics_inp must be set
*.vcham file structure

In this file all parameters for the vibronic-coupling Hamiltonian are written down. All this parameters are necassary to describe the potential energy surfaces. The unit is eV.
The file is structurated into different parts. On the beginning all parameters are list: frequencies, equilibrium point energies at Q0, linear on- and off-diagonal constants, bilinear on- and off-diagonal constants.
Then follows a section where the energyminimas are listed.
At the end you can get the information which parameters were calculated. 0 means not and 1 means yes.

*.par file structure

This file includes the same parameters as in the *.vcham file. But this file is in a format, that can be used as MCTDH input.

*.op file structure

This file includes the same parameters as in the *.par file. This file aslo includes the hamiltonian section, in a format that can be used as MCTDH input.

*.log file structure

Some useful information are given in the *.log file.

Example: Cr(CO)5
As an example, crco5_vchfit0.inp uses the information in the crco5_trans0.info database (see VCTRANS documentation) to get the parameters. No fit is made, but the derivative coupling and gradient difference information from CASSCF calculations are used to obtain the parameters. The results are in the file crco5_vchfit0.vcham and crco5_vchfit0.log

This calculation is then used as the initial guess in crco5_vchfit1.inp a fit of the model to the data in the full database crco5_trans1.info . The results are written to the files crco5_vchfit1.vcham and crco5_vchfit1.log

The next step is to improve the fit. crco5_vchfit2.inp uses crco5_vchfit1.vcham as an initial guess, but includes an energyweight to improve the fit in the minima and also optimises the on-diagonal 2nd-order "gamma" constants. Results are written to the files crco5_vchfit2.vcham and crco5_vchfit2.log

An improved fit is then obtained using crco5_vchfit3.inp. This uses crco5_vchfit2.vcham as an initial guess, and also includes quartic on-diagonal terms. Results are written to the files crco5_vchfit3.vcham and crco5_vchfit3.log As the keyword mctdh_op had been used, the MCTDH operator file crco5_vchfit3.op has also been written.