For interpretive purposes, one might want to look at the wavefunctions produced by the xsim package. If xsim is used for a single electronic state, these are vibrational wavefunctions; if it is used by more than one state, these are vibronic wavefunctions associated with a Born-Huang expansion. Before discussing what these "mean" and how to look at them, here are the nuts and bolts of generating them.
1. Obtain a vibronic wavefunction, either by doing a Davidson diagonalization or by appropriate linear combinations of Lanczos trial vectors (see below). If the wavefunction is generated by the Davidson procedure, it will be found in the file called GUESSVECTOR. Rename this file to plotfile.in, and go to (3) below.
2. It is often hard/impossible to converge excited state roots to get precise eigenvectors. However, there is a provision to generate (quite good) approximate wavefunctions by expanding the (near) eigenvector in the basis of the tridiagonal Lanczos representation. In order to do this:
a) You must store the Lanczos trial vectors and the eigenvectors of the ultimate tridiagonal matrix. To do this, simply add the keyword
to your xsim input file. However: BE CAREFUL. A typical xsim calculation for me will include perhaps 10,000,000 basis functions and 1,000 Lanczos iterations. Such a calculation would result in the storage of 1000 vectors, which are all 80 Mb in length, for a total of 80 Gb (the additional storage of the 1000 eigenvectors of the 1000x1000 tridiagonal matrix is a trivial addition here). Using Storevectors with a large basis set and many Lanczos iterations is perhaps a good way to get you on the wrong side of your friendly system manager. This information will be stored in the LEVECS file (the eigenvectors of the tridiagonal matrix and the associated eigenvalues) and LBASIS, which contains the Lanczos trial vectors. Now, the real eigenvectors of the Hamiltonian in the full direct-product Born-Huang basis are approximated by linear combinations of the Lanczos trial vectors, with the weights specified by the appropriate eigenvector.
b) The easiest way to construct the Lanczos approximation to the full eigenvectors is to add the
keyword to your xsim input file. This record will look something like this:
-0.395698072715971 1.372600825906208E-002 558.840893830490
-0.336615182517430 4.281480354483751E-002 1035.37394472682
-0.307337923582097 2.572448753157076E-002 1271.50967666976
-0.298430595114833 2.865920756915017E-002 1343.35173442247
-0.273975449003787 3.160906765614277E-002 1540.59471538112
-0.256027361967664 1.512601468720900E-003 1685.35501137096
the entry immediately after Plotroots means that six wavefunctions will be constructed. The following lines (which are taken from the fort.20 file produced by xsim) give the eigenvalues (eV in left-most column), intensities (in central column) and offset from the lowest root (right-hand column, in cm-1). The xsimconst program uses this information to find the roots in the LBASIS and LEVECS files that correspond to these states, and then proceeds to construct the remaining vector.
In general, roots that have vanishing (or nearly vanishing intensity) converge very slowly in the Lanczos algorithm. Constructing "wavefunctions" for such states by the method described here is not reliable and can, in fact, produce what might charitably be called "garbage". If you want to plot such a state, it is wise choose a Lanczos seed vector (through the Seed keyword in the xsim input file) which will have a significant project upon the state of interest.
c) With the Plotroots input shown above, the actual Lanczos approximations to the Hamiltonian eigenvectors are generated by
xsimconst < input
where input is again the xsim input file that contains the Plotroots record.
d) After running xsimconst, the directory willl contain files with names like:
outfile00558.84 outfile01271.50 outfile01540.59
outfile01035.37 outfile01343.35 outfile01685.35
which are the wavefunctions for the states with the associated energy. These can now be rennamed to plotfile.in (as above), and you can now proceed to (3) below.
3. Generate a plot with the xvibplot executable:
xvibplot < input
where input is the usual xsim input file that was used to generate the wavefunction. In order to make a plot, one additional keyword must be included:
which means that the wavefunction will be projected onto the x,y plane where x and y are the numbers of the corresponding xsim input file coordinates. For the moment (although this will be changed to be optional in the future), all other coordinates are set to 0. For example, if you have a six mode system, and your Plotcoord record is
then running xvibplot will generate wavefunctions over a grid of points (0 x 0 0 0 z 0). By default, the grid over x and z extends from -10 to 10 in each coordinate, with sampling every 0.4 (dimensionless normal coordinates are used throughout by xsim).
After xvibplot runs, a number of files will be found. The files psia, psib, psic ... are the vibrational parts of the wavefunction associated with the first, second, third, etc. electronic state in the basis (which will, in general have different symmetries, as the overall symmetry of the vibronic wavefunction is the product of the vibrational symmetry and the electronic symmetries, all of which will be the same throughout the Born-Huang expansion). The file psisquared is exactly what is implied by the name; it is the probability density of the vibronic state in the nuclear coordinates.
All of these functions can be nicely plotted with gnuplot as follows:
splot 'psia' w points palette
which will show a perspective view. It is perhaps most easy to "interpret" the wavefunctions from an "overhead" view, which can be generated by adding the "set view" gnuplot directive:
set view 0,0
splot 'psia' w points palette
Of course, the same recipe works for psib, psic ... and psisquared. The last of these is generally the most useful one.