using-vcham

How to use VCHam - a program to calculate the vibronic coupling Hamiltonian.

VCHam toolbox

vcpnt101, generates quantum chemistry input files (e.g. g09) for point along a chosen normal mode

vctrans101, reads data from quantum chemistry output files (generated by vcpnt) and formats for use in vchfit

vchfit101, fits curves to potential energy curves along a normal mode

Example usage

In this example there are 30 normal modes and 5 electronic states.

vcpnt

vcpnt101 pnts.inp

will generate input files along a chosen normal mode according to pnts.inp; Here is an example pnts.inp file:

file0 = structure1.log
nmodes = 30
fileout = structure1
opening = open.dat
closing = close.dat
formout = gaussian
along = 01 10 1.0

end-input

Where,

• file0, the output for a R0 frequency calculation (e.g. opt, freq in g09);

• nmodes, the total number of normal modes (3N-6 for non-linear molecules)

• fileout, the title of the files to be generated by vcpnt. They will be this title followed by the mode number and the number of deviations left or right from the R0 position (as read from file0 = structure1.log) separated by underscores.

• opening, this specifies a special file containing the opening for a g09 template input file (before the geometry section).

• closing, this specifies a special file containing the closing for a g09 template input file (after the geometry section).

• formout, specifies the format of the files generated by vcpnt; in this case gaussian

• along, specifies the mode number, number of deviations (left and right) from the R0 position, and the spacing (not sure of units, normal mode coordinates?) 1.0 is a standard value

Then, run the generated g09 calculations (will take a while).

vctrans

Next, use vctrans to transform the data from the g09 output files to nice format and for use in vchfit,

vctrans101 trans.inp

As this must be done for every mode, use a template file such as 'trans_xx.inp', which contains 'xx' in place of where the mode number should be. Here's an example 'trans_xx.inp':

file0 = hpfreq.log
nmodes = 156
orient=Standard
abinitiotype =  tddft            #ttddft for triplets
energy0 = -1288.75520321
states = 1,2,3,4,5
order = 0

files
pnts_xx_000.log 1
pnts_xx_000.log 2
pnts_xx_000.log 3
pnts_xx_000.log 4
pnts_xx_000.log 5
end-files

datasets
 qxx.set
end-datasets

end-input

Then we can just,

for i in {1..30}
do
sed "s/xx/$i/g" trans_xx.inp > trans_$i.inp
done

The trans.inp files reference a 'q.set' file which is a list of g09 log files which VCHam reads from. Similarly, just use a template file 'qxx.set', e.g:

pnts_xx_10l.log 1
pnts_xx_10l.log 2
pnts_xx_10l.log 3
pnts_xx_10l.log 4
pnts_xx_10l.log 5
pnts_xx_09l.log 1
pnts_xx_09l.log 2
pnts_xx_09l.log 3
pnts_xx_09l.log 4
pnts_xx_09l.log 5
...
pnts_xx_01l.log 1
pnts_xx_01l.log 2
pnts_xx_01l.log 3
pnts_xx_01l.log 4
pnts_xx_01l.log 5
pnts_xx_10r.log 1
pnts_xx_10r.log 2
pnts_xx_10r.log 3
pnts_xx_10r.log 4
pnts_xx_10r.log 5
pnts_xx_09r.log 1
pnts_xx_09r.log 2
pnts_xx_09r.log 3
pnts_xx_09r.log 4
pnts_xx_09r.log 5
...
pnts_xx_01r.log 1
pnts_xx_01r.log 2
pnts_xx_01r.log 3
pnts_xx_01r.log 4
pnts_xx_01r.log 5
end-files

Generate q1.set, q2.set, etc.,

for i in {1..30}
do
sed "s/xx/$i/g" qxx.set > q$i.set
done

Then run vctrans on all the trans.inp files,

for i in trans*.inp
do
vctrans101 $i
done

To generate a 'trans.info' file for each mode. They can be plotted with vcplot, e.g.

vcplot101 trans_1.inp

vcfit

To fit the PECs along each normal mode, Fit-x-1.inp, Fit-x-2.inp and Fit-x-3.inp files must be used. Here is an example Fit-x-2.inp:

infofile = trans_x.info
adiab2_mod
mctdh_op

iter = 64

fit_select = x

symmode = A

cpmode = 1 2 B  cpmode = 1 3 A  cpmode = 1 4 B  cpmode = 1 5 B
cpmode = 2 3 B  cpmode = 2 4 A  cpmode = 2 5 A
cpmode = 3 4 B  cpmode = 3 5 B
cpmode = 4 5 A

diabatic_function
end-define

kappafit
lambdafit
gammafit

guess=fitx-1.vcham
#use_cons

constraints-section
end-constraints-section

end-input

This file involves the cpmode section which is the coupling between states. In this case, we have a molecule with C2 symmetry. The ground state (1) is always of symmetry A; in this case, the 1st excited singlet is of irrep B so the overall irrep is B because AuB = B (and AuA = A, and BuB = A). Thus, we have 'cpmode = 1 2 B', that is, irrep B between states 1 and 2 (ground and 1st excited singlet). Note the cpmode section is there only to speed up the fitting as it allows certain matrix elements to be 0 due to symmetry. We can omit the entire section and proceed, albeit slower computationally. A molecule of C1 symmetry for example has no advantage of cpmode. If that were the case, just omit the cpmode section.

We actually use three of these for different fitting algorithms, eventually obtaining a good final fit. They differ by lambdafit (fit01), then gammafit (fit02), then we fill out the diabatic_function section for mode x to use quartic terms (fit03),

diabatic_function
x 1 quartic
x 2 quartic
x 3 quartic
x 4 quartic
x 5 quartic
end-define

Also 'guess=' becomes uncommented for fit02 and fit03, using the corresponding previous fit as its guess. To fit we have to generate fit01, fit02, and fit03 for each mode e.g. fit-1-1.inp, fit-1-2.inp, fit-1-3.inp, fit-2-1.inp, fit-2-2.inp, fit-2-3.inp, ... , up to fit-30-1.inp, fit-30-2.inp, fit-30-3.inp. Use a template file and sed again,

for i in {1..30}
do
sed "s/xx/$i/g" fit_xx-1.inp > fit_$i-1.inp
sed "s/xx/$i/g" fit_xx-2.inp > fit_$i-2.inp
sed "s/xx/$i/g" fit_xx-3.inp > fit_$i-3.inp
done

Then run the fitting e.g.,

for i in fit*.inp 
do
vchfit101 $i
done

This works because fit_xx-1.inp is listed before fit_xx-2.inp which are both before fit_xx-3.inp in bash; try it by

echo fit*.inp

The fittings must be in order 1,2,3, because fit2 uses info from fit1 as a first guess, and fit3 uses fit2 as a first guess.

Finally, they can be plotted again including the fitted lines. Usually we plot fit3, the most accurate fit,

vcplot101 fit-1-3