How to use VCHam - a program to calculate the vibronic coupling Hamiltonian.
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
In this example there are 30 normal modes and 5 electronic states.
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).
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
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