In C_{2v}, however, the functions are distributed into the four representations of the group and therefore different symmetry representations can be mixed. The next table lists the distribution of the functions in C_{2v} and the symmetry of the corresponding orbitals in .
In symmetry we find both and orbitals. When the calculation is performed in C_{2v} symmetry all the orbitals of symmetry can mix because they belong to the same representation, but this is not correct for . The total symmetry must be kept and therefore the orbitals should not be allowed to rotate and mix with the orbitals. The same is true in the and symmetries with the and orbitals, while in symmetry this problem does not exist because it has only orbitals (with a basis set up to f functions). The tool to restrict possible orbital rotations is the option SUPSym in the RASSCF program. It is important to start with clean orbitals belonging to the actual symmetry, that is, without unwanted mixing. But the problems with the symmetry are not solved with the SUPSym option only. Orbitals belonging to different components of a degenerate representation should also be equivalent. For example: the orbitals in and symmetries should have the same shape, and the same is true for the orbitals in and symmetries. This can only be partly achieved in the RASSCF code. The input option AVERage will average the density matrices for representations and ( and orbitals), thus producing equivalent orbitals. The present version does not, however, average the orbital densities in representations and (note that this problem does not occur for electronic states with an equal occupation of the two components of a degenerate set, for example states). A safe way to obtain totally symmetric orbitals is to reduce the symmetry to C_{1} (or C_{s} in the homonuclear case) and perform a stateaverage calculation for the degenerate components. We need an equivalence table to know the correspondence of the symbols for the functions in MOLCAS to the spherical harmonics (SH):
We begin by performing a SCF calculation and analyzing the resulting orbitals. The employed bond distance is close to the experimental equilibrium bond length for the ground state [209]. Observe in the following SEWARD input that the symmetry generators, planes yz and xz, lead to a C_{2v} representation. In the SCF input we have used the option OCCNumbers which allows specification of occupation numbers other than 0 or 2. It is still the closed shell SCF energy functional which is optimized, so the obtained SCF energy has no physical meaning. However, the computed orbitals are somewhat better for open shell cases as NiH. The energy of the virtual orbitals is set to zero due to the use of the IVO option. The order of the orbitals may change in different computers and versions of the code.
In difficult situations it can be useful to employ the AUFBau option of the SCF program. Including this option, the subsequent classification of the orbitals in the different symmetry representations can be avoided. The program will look for the lowestenergy solution and will provide with a final occupation. This option must be used with caution. It is only expected to work in clear closedshell situations. We have only printed the orbitals most relevant to the following discussion. Starting with symmetry 1 () we observe that the orbitals are not mixed at all. Using a basis set contracted to Ni 5s4p3d1f / H 3s2p in symmetry we obtain 18 molecular orbitals (combinations from eight atomic s functions, six p_{z} functions, three d_{z2} functions, and one f_{z3} function) and four orbitals (from three d_{x2y2} functions and one f_{z(x2y2)} function). Orbitals 6, 10, 13, and 18 are formed by contributions from the three d_{x2y2} and one f_{z(x2y2)} functions, while the contributions of the remaining harmonics are zero. These orbitals are orbitals and should not mix with the remaining orbitals. The same situation occurs in symmetries and (2 and 3) but in this case we observe an important mixing among the orbitals. Orbitals 7 and 7 have main contributions from the harmonics 4f3+ (f_{x3}) and 4f3 (f_{y3}), respectively. They should be pure orbitals and not mix at all with the remaining orbitals. The first step is to evaluate the importance of the mixings for future calculations. Strictly, any kind of mixing should be avoided. If g functions are used, for instance, new contaminations show up. But, undoubtedly, not all mixings are going to be equally important. If the rotations occur among occupied or active orbitals the influence on the results is going to be larger than if they are high secondary orbitals. NiH is one of these cases. The ground state of the molecule is . It has two components and we can therefore compute it by placing the single electron in the d_{xy} orbital (leading to a state of symmetry in C_{2v}) or in the d_{x2y2} orbital of the symmetry. Both are orbitals and the resulting states will have the same energy provided that no mixing happens. In the symmetry no mixing is possible because it is only composed of orbitals but in symmetry the and orbitals can rotate. It is clear that this type of mixing will be more important for the calculation than the mixing of and orbitals. However it might be necessary to prevent it. Because in the SCF calculation no high symmetry restriction was imposed on the orbitals, orbitals 2 and 4 of the and symmetries have erroneous contributions of the 4f3+ and 4f3 harmonics, and they are occupied or active orbitals in the following CASSCF calculation. To use the supersymmetry (SUPSym) option we must start with proper orbitals. In this case the orbitals are symmetry adapted (within the printed accuracy) but not the and orbitals. Orbitals 7 and 7 must have zero coefficients for all the harmonics except for 4f3+ and 4f3, respectively. The remaining orbitals of these symmetries (even those not shown) must have zero in the coefficients corresponding to 4f3+ or 4f3. To clean the orbitals the option CLEAnup of the RASSCF program can be used. Once the orbitals are properly symmetrized we can perform CASSCF calculations on different electronic states. Deriving the types of the molecular electronic states resulting from the electron configurations is not simple in many cases. In general, for a given electronic configuration several electronic states of the molecule will result. Wigner and Witmer derived rules for determining what types of molecular states result from given states of the separated atoms. In chapter VI of reference [211] it is possible to find the tables of the resulting electronic states once the different couplings and the Pauli principle have been applied. In the present CASSCF calculation we have chosen the active space (3d, 4d, , ) with all the 11 valence electrons active. If we consider 4d and as weakly occupied correlating orbitals, we are left with 3d and (six orbitals), which are to be occupied with 11 electrons. Since the bonding orbital (composed mainly of Ni 4s and H 1s) will be doubly occupied in all low lying electronic states, we are left with nine electrons to occupy the 3d orbitals. There is thus one hole, and the possible electronic states are: , , and , depending on the orbital where the hole is located. Taking Table into account we observe that we have two lowlying electronic states in symmetry 1 (A_{1}): and , and one in each of the other three symmetries: in symmetries 2 (B_{1}) and 3 (B_{2}), and in symmetry 4 (A_{2}). It is not immediately obvious which of these states is the ground state as they are close in energy. It may therefore be necessary to study all of them. It has been found at different levels of theory that the NiH has a ground state [209]. We continue by computing the ground state. The previous SCF orbitals will be the initial orbitals for the CASSCF calculation. First we need to know in which C_{2v} symmetry or symmetries we can compute a state. In the symmetry tables it is determined how the species of the linear molecules are resolved into those of lower symmetry (depends also on the orientation of the molecule). In Table is listed the assignment of the different symmetries for the molecule placed on the z axis. The state has two degenerate components in symmetries and . Two CASSCF calculations can be performed, one computing the first root of symmetry and the second for the first root of symmetry. The RASSCF input for the state of symmetry would be:
The corresponding input for symmetry will be identical except
for the SYMMetry keyword
In the RASSCF inputs the CLEAnup option will take the initial orbitals (SCF here) and will place zeroes in all the coefficients of orbitals 6, 10, 13, and 18 in symmetry 1, except in coefficients 13, 14, 15, and 17. Likewise all coefficients 13, 14, 15, and 17 of the remaining orbitals will be set to zero. The same procedure is used in symmetries and . Once cleaned, and because of the SUPSymmetry option, the orbitals 6, 10, 13, and 18 of symmetry will only rotate among themselves and they will not mix with the remaining orbitals. The same holds true for orbitals 7 and 7 in their respective symmetries. Orbitals can change order during the calculation. MOLCAS incorporates a procedure to check the nature of the orbitals in each iteration. Therefore the right behavior of the SUPSym option is guaranteed during the calculation. The procedure can have problems if the initial orbitals are not symmetrized properly. Therefore, the output with the final results should be checked to compare the final order of the orbitals and the final labeling of the SUPSym matrix. The AVERage option would average the density matrices of symmetries 2 and 3, corresponding to the and symmetries in . In this case it is not necessary to use the option because the two components of the degenerate sets in symmetries and have the same occupation and therefore they will have the same shape. The use of the option in a situation like this ( and states) leads to convergence problems. The symmetry of the orbitals in symmetries 2 and 3 is retained even if the AVERage option is not used. The output for the calculation on symmetry 4 () contains the following lines:
The state is mainly (weight 96%) described by a single configuration (configuration number 15834) which placed one electron on the first active orbital of symmetry 4 () and the remaining electrons are paired. A close look to this orbital indicates that is has a coefficient .9989 in the first 3d2 (3d_{xy}) function and small coefficients in the other functions. This results clearly indicate that we have computed the state as the lowest root of that symmetry. The remaining configurations have negligible contributions. If the orbitals are properly symmetrized, all configurations will be compatible with a electronic state. The calculation of the first root of symmetry 1 () results:
We obtain the same energy as in the previous calculation. Here the dominant configuration places one electron on the first active orbital of symmetry 1 (). It is important to remember that the orbitals are not ordered by energies or occupations into the active space. This orbital has also the coefficient .9989 in the first 3d2 (3d_{x2y2}) function. We have then computed the other component of the state. As the orbitals in different C_{2v} symmetries are not averaged by the program it could happen (not in the present case) that the two energies differ slightly from each other. The consequences of not using the SUPSym option are not extremely severe in the present example. If you perform a calculation without the option, the obtained energy is:
As it is a broken symmetry solution the energy is lower than in the other case. This is a typical behavior. If we were using an exact wave function it would have the right symmetry properties, but approximated wave functions do not necessarily fulfill this condition. So, more flexibility leads to lower energy solutions which have broken the orbital symmetry. If in addition to the state we want to compute the lowest state we can use the adapted orbitals from any of the state calculations and use the previous RASSCF input without the CLEAnup option. The orbitals have not changed place in this example. If they do, one has to change the labels in the SUPSym option. The simplest way to compute the lowest excited state is having the unpaired electron in one of the orbitals because none of the other configurations, or , leads to the term. However, there are more possibilities such as the configuration ; three nonequivalent electrons in three orbitals. In actuality the lowest state must be computed as a doublet state in symmetry A_{1}. Therefore, we set the symmetry in the RASSCF to 1 and compute the second root of the symmetry (the first was the state):
Of course the SUPSym option must be maintained. The use of CIROot indicates that we are computing the second root of that symmetry. The obtained result:
As we have used two as the dimension of the CI matrix employed in the CI Davidson procedure we obtain the wave function of two roots, although the optimized root is the second. Root 1 places one electron in the first active orbital of symmetry one, which is a 3d2+ (3d_{x2y2}) orbital. Root 2 places the electron in the second active orbital, which is a orbital with a large coefficient (.9639) in the first 3d0 (3d_{z2}) function of the nickel atom. We have therefore computed the lowest state. The two states resulting from the configuration with the three unpaired electrons is higher in energy at the CASSCF level. If the second root of symmetry had not been a state we would have to study higher roots of the same symmetry. It is important to remember that the active orbitals are not ordered at all within the active space. Therefore, their order might vary from calculation to calculation and, in addition, no conclusions about the orbital energy, occupation or any other information can be obtained from the order of the active orbitals. We can compute also the lowest excited state. The simplest possibility is having the configuration , which only leads to one state. The unpaired electron will be placed in either one or one orbital. That means that the state has two degenerate components and we can compute it equally in both symmetries. There are more possibilities, such as the configuration or the configuration . The resulting state will always have two degenerate components in symmetries and , and therefore it is the wave function analysis which gives us the information of which configuration leads to the lowest state. For NiH it turns out to be non trivial to compute the state. Taking as initial orbitals the previous SCF orbitals and using any type of restriction such as the CLEAnup, SUPSym or AVERage options lead to severe convergence problems like these:
The calculation, however, converges in an straightforward way if none of those tools are used:
The (and ) orbitals, both in symmetries and , are, however, differently occupied and therefore are not equal as they should be:
Therefore what we have is a symmetry broken solution. To obtain a solution which is not of broken nature the and orbitals must be equivalent. The tool to obtain equivalent orbitals is the AVERage option, which averages the density matrices of symmetries and . But starting with any of the preceding orbitals and using the AVERage option lead again to convergence problems. It is necessary to use better initial orbitals; orbitals which have already equal orbitals in symmetries and . One possibility is to perform a SCF calculation on the NiH cation explicitly indicating occupation one in the two higher occupied orbitals (symmetries 2 and 3):
It can take some successive steps to obtain a converged calculation using the CLEAnup, SUPSym, and AVERage options. The calculation with a single root did not converge clearly. We obtained, however, a converged result for the lowest state of NiH by computing two averaged CASSCF roots and setting a weight of 90% for the first root using the keyword:
The energy of the different states (only the first one shown above) is printed on the top of their configuration list. The converged energy is simply an average energy. The occupation numbers obtained in the section of the RASSCF output printed above are the occupation numbers of the natural orbitals of the corresponding root. They differ from the occupation numbers printed in the molecular orbital section where we have pseudonatural molecular orbitals and average occupation numbers. On top of each of the valence orbitals an average occupation close to 1.5e will be printed; this is a consequence of the the averaging procedure. The results obtained are only at the CASSCF level. Additional effects have to be considered and included. The most important of them is the dynamical correlation effect which can be added by computing, for instance, the CASPT2 energies. The reader can find a detailed explanation of the different approaches in ref. [209], and a careful discussion of their consequences and solutions in ref. [212]. We are going, however, to point out some details. In the first place the basis set must include up to g functions for the transition metal atom and up to d functions for the hydrogen. Relativistic effects must be taken into account, at least in a simple way as a first order correction. The keyword RELInt must be then included in the SEWARD input to compute the massvelocity and oneelectron Darwin contact term integrals and obtain a firstorder correction to the energy with respect to relativistic effects at the CASSCF level in the RASSCF output. Scalar relativistic effects can be also included according the DouglasKroll or the BaryszSadlejSnijders transformations, as it will be explained in section . The CASPT2 input needed to compute the secondorder correction to the energy will include the number of the CASSCF root to compute. For instance, for the first root of each symmetry:
The number of frozen orbitals taken by CASPT2 will be that specified in the RASSCF input except if this is changed in the CASPT2 input. In the perturbative step we have frozen all the occupied orbitals except the active ones. This is motivated by the desire to include exclusively the dynamical correlation related to the valence electrons. In this way we neglect correlation between core electrons, named corecore correlation, and between core and valence electrons, named corevalence correlation. This is not because the calculation is smaller but because of the inclusion of those type of correlation in a calculation designed to treat valence correlation is an inadequate approach. Corecore and corevalence correlation requires additional basis functions of the same spatial extent as the occupied orbitals being correlated, but with additional radial and angular nodes. Since the spatial extent of the core molecular orbitals is small, the exponents of these correlating functions must be much larger than those of the valence optimized basis sets. The consequence is that we must avoid the inclusion of the core electrons in the treatment in the first step. Afterwards, the amount of correlation introduced by the core electrons can be estimated in separated calculations for the different states and those effects added to the results with the valence electrons. Corevalence correlation effects of the 3s and 3p nickel shells can be studied by increasing the basis set flexibility by uncontracting the basis set in the appropriate region. There are different possibilities. Here we show the increase of the basis set by four s, four p, and four d functions. f functions contribute less to the description of the 3s and 3p shells and can be excluded. The uncontracted exponents should correspond to the region where the 3s and 3p shells present their density maximum. Therefore, first we compute the absolute maxima of the radial distribution of the involved orbitals, then we determine the primitive gaussian functions which have their maxima in the same region as the orbitals and therefore which exponents should be uncontracted. The final basis set will be the valence basis set used before plus the new added functions. In the present example the SEWARD input can be:
We have used a special format to include the additional functions. We include the additional 4s4p4d functions for the nickel atom. The additional basis set input must use a dummy label (Nix here), the same coordinates of the original atom, and specify a CHARge equal to zero, whether in an Inline basis set input as here or by specifically using keyword CHARge. It is not necessary to include the basis set with the Inline format. A library can be created for this purpose. In this case the label for the additional functions could be:
and a proper link to AUXLIB should be included in the script (or in the input if one uses AUTO). Now the CASPT2 is going to be different to include also the correlation related to the 3s,3p shell of the nickel atom. Therefore, we only freeze the 1s,2s,2p shells:
A final effect one should study is the basis set superposition error (BSSE). In many cases it is a minor effect but it is an everpresent phenomenon which should be investigated when high accuracy is required, especially in determining bond energies, and not only in cases with weakly interacting systems, as is frequently believed. The most common approach to estimate this effect is the counterpoise correction: the separated fragment energies are computed in the total basis set of the system. For a discussion of this issue see Refs. [212,213]. In the present example we would compute the energy of the isolated nickel atom using a SEWARD input including the full nickel basis set plus the hydrogen basis set in the hydrogen position but with the charge set to zero. And then the opposite should be done to compute the energy of isolated hydrogen. The BSSE depends on the separation of the fragments and must be estimated at any computed geometry. For instance, the SEWARD input necessary to compute the isolated hydrogen atom at a given distance from the ghost nickel basis set including core uncontracted functions is:
The results obtained at the CASPT2 level are close to those obtained by MRCI+Q and ACPF treatments but more accurate. They match well with experiment. The difference is that all the configuration functions (CSFs) of the active space can be included in CASPT2 in the zerothorder references for the secondorder perturbation calculation [209], while the other methods have to restrict the number of configurations. Calculations of linear molecules become more and more complicated when the number of unpaired electrons increases. In the following sections we will discuss the more complicated situation occurring in the Ni_{2} molecule.

Symmetry  Spherical harmonics  
s  d_{z2}  
p_{z}  f_{z3}  
d_{xz}  d_{yz}  
p_{x}  p_{y}  f_{x(z2y2)}  f_{y(z2x2)}  
d_{x2y2}  d_{xy}  
f_{xyz}  f_{z(x2y2)}  
f_{x3}  f_{y3}  
^{a}Functions placed on the symmetry center. 
Table classifies the functions and orbitals into the symmetry representations of the D_{2h} symmetry. Note that in table subindex b stands for bonding combination and a for antibonding combination.
Symm.^{b}  Spherical harmonics (orbitals in )  
(1)  s_{b} ()  p_{zb} ()  d_{z2b} ()  d_{x2y2b} ()  f_{z3b} ()  f_{z(x2y2)b} () 
(2)  p_{xb} ()  d_{xzb} ()  f_{x(z2y2)b} ()  f_{x3b} ()  
(3)  p_{yb} ()  d_{yzb} ()  f_{y(z2x2)b} ()  f_{y3b} ()  
(4)  d_{xyb} ()  f_{xyzb} ()  
(5)  s_{a} ()  p_{za} ()  d_{z2a} ()  d_{x2y2a} ()  f_{z3a} ()  f_{z(x2y2)a} () 
(6)  p_{ya} ()  d_{yza} ()  f_{y(z2x2)a} ()  f_{y3a} ()  
(7)  p_{xa} ()  d_{xza} ()  f_{x(z2y2)a} ()  f_{x3a} ()  
(8)  d_{xya} ()  f_{xyza} ()  
^{a}Subscripts a and b refer to the bonding and antibonding combination of the AO's, respectively.  
^{b}In parenthesis the number of the symmetry in MOLCAS. Note that the number and order of the  
symmetries depend on the generators and the orientation of the molecule. 
The order of the symmetries, and therefore the number they have in MOLCAS, depends on the generators used in the SEWARD input. This must be carefully checked at the beginning of any calculation. In addition, the orientation of the molecule on the cartesian axis can change the labels of the symmetries. In Table for instance we have used the order and numbering of a calculation performed with the three symmetry planes of the D_{2h} point group (X Y Z in the SEWARD input) and the z axis as the intermolecular axis (that is, x and y are equivalent in D_{2h}). Any change in the orientation of the molecule will affect the labels of the orbitals and states. In this case the orbitals will belong to the , , , and symmetries. For instance, with x as the intermolecular axis and will be replaced by and , respectively, and finally with y as the intermolecular axis , , , and would be the orbitals.
It is important to remember that MOLCAS works with symmetry adapted basis functions. Only the symmetry independent atoms are required in the SEWARD input. The remaining ones will be generated by the symmetry operators. This is also the case for the molecular orbitals. MOLCAS will only print the coefficients of the symmetry adapted basis functions.
The necessary information to obtain the complete set of orbitals
is contained in the SEWARD output. Consider the case of the symmetry:
**************************************************
******** Symmetry adapted Basis Functions ********
**************************************************
~
Irreducible representation : ag
Basis function(s) of irrep:
~
Basis Label Type Center Phase Center Phase
1 C 1s0 1 1 2 1
2 C 1s0 1 1 2 1
3 C 1s0 1 1 2 1
4 C 1s0 1 1 2 1
5 C 2pz 1 1 2 1
6 C 2pz 1 1 2 1
7 C 2pz 1 1 2 1
8 C 3d0 1 1 2 1
9 C 3d0 1 1 2 1
10 C 3d2+ 1 1 2 1
11 C 3d2+ 1 1 2 1
12 C 4f0 1 1 2 1
13 C 4f2+ 1 1 2 1
The previous output indicates that symmetry adapted basis function 1, belonging to the representation, is formed by the symmetric combination of a s type function centered on atom C and another s type function centered on the redundant center 2, the second carbon atom. Combination s+s constitutes a bonding type orbital. For the p_{z} function however the combination must be antisymmetric. It is the only way to make the p_{z} orbitals overlap and form a bonding orbital of symmetry. Similar combinations are obtained for the remaining basis sets of the and other symmetries.
The molecular orbitals will be combinations of these symmetry adapted functions. Consider the orbitals:
SCF orbitals
~
Molecular orbitals for symmetry species 1
~
ORBITAL 1 2 3 4 5 6
ENERGY 11.3932 1.0151 .1138 .1546 .2278 .2869
OCC. NO. 2.0000 2.0000 .0098 .0000 .0000 .0000
~
1 C 1s0 1.4139 .0666 .0696 .2599 .0626 .0000
2 C 1s0 .0003 1.1076 .6517 1.0224 .4459 .0000
3 C 1s0 .0002 .0880 .2817 .9514 .0664 .0000
4 C 1s0 .0000 .0135 .0655 .3448 .0388 .0000
5 C 2pz .0006 .2581 1.2543 1.1836 .8186 .0000
6 C 2pz .0000 .1345 .0257 2.5126 1.8556 .0000
7 C 2pz .0005 .0192 .0240 .7025 .6639 .0000
8 C 3d0 .0003 .0220 .0005 .9719 .2430 .0000
9 C 3d0 .0001 .0382 .0323 .8577 .2345 .0000
10 C 3d2+ .0000 .0000 .0000 .0000 .0000 .7849
11 C 3d2+ .0000 .0000 .0000 .0000 .0000 .7428
12 C 4f0 .0002 .0103 .0165 .0743 .0081 .0000
13 C 4f2+ .0000 .0000 .0000 .0000 .0000 .0181
In MOLCAS outputs only 13 coefficients for orbital are going to be printed because they are the coefficients of the symmetry adapted basis functions. If the orbitals were not composed by symmetry adapted basis functions they would have, in this case, 26 coefficients, two for type of function (following the scheme observed above in the SEWARD output), symmetrically combined the s and d functions and antisymmetrically combined the p and f functions.
To compute electronic states using the D_{2h} symmetry we need to go to the symmetry tables and determine how the species of the linear molecules are resolved into those of lower symmetry (this depends also on the orientation of the molecule [211]). Table lists the case of a linear molecule with z as the intermolecular axis.
State symmetry  State symmetry D_{2h}  
A_{g}  
B_{1u}  
B_{1g}  
A_{u}  
B_{2g} + B_{3g}  
B_{2u} + B_{3u}  
A_{g} + B_{1g}  
A_{u} + B_{1u}  
B_{2g} + B_{3g}  
B_{2u} + B_{3u}  
A_{g} + B_{1g}  
A_{u} + B_{1u} 
To compute the ground state of C_{2}, a state, we will
compute a singlet state of symmetry A_{g} (1 in this context).
The input files for a CASSCF calculation on the C_{2} ground state
will be:
&SEWARD &END
Title
C2
Symmetry
X Y Z
Basis set
C.ANOL...4s3p2d1f.
C .00000000 .00000000 1.4
End of basis
End of input
&SCF &END
Title
C2
ITERATIONS
40
Occupied
2 1 1 0 2 0 0 0
End of input
&RASSCF &END
Title
C2
Nactel
4 0 0
Spin
1
Symmetry
1
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
*Average
*2 2 3 6 7
Supsymmetry
1
3 6 9 11
1
1 6
1
1 6
0
1
3 5 8 12
1
1 6
1
1 6
0
Iter
50,25
Lumorb
End of input
In this case the SCF orbitals are already clean symmetry adapted orbitals (within the printed accuracy). We can then directly use the SUPSym option. In symmetries and we restrict the rotations among the and the orbitals, and in symmetries , , , and the rotations among and orbitals. Additionally, symmetries and and symmetries and are averaged, respectively, by using the AVERage option. They belong to the and representations in , respectively.
A detailed explanation on different CASSCF calculations on the C_{2} molecule and their states can be found elsewhere [73]. Instead we include here an example of how to combine the use of UNIX shell script commands with MOLCAS as a powerful tool.
The following example computes the transition dipole moment for the transition from the state to the state in the C_{2} molecule. This transition is known as the Phillips bands [211]. This is not a serious attempt to compute this property accurately, but serves as an example of how to set up an automatic calculation. The potential curves are computed using CASSCF wavefunctions along with the transition dipole moment.
Starting orbitals are generated by computing a CI wavefunction once and using the natural orbitals. We loop over a set of distances, compute the CASSCF wave functions for both states and use RASSI to compute the TDMs. Several UNIX commands are used to manipulate input and output files, such as grep, sed, and the awk language. For instance, an explicit 'sed' is used to insert the geometry into the seward input; the final CASSCF energy is extracted with an explicit 'grep', and the TDM is extracted from the RASSI output using an awk script. We are not going to include the awk scripts here. Other tools can be used to obtain and collect the data.
In the first script, when the loop over geometries is done, four files are available: geom.list (contains the distances), tdm.list (contains the TDMs), e1.list (contains the energy for the state), and e2.list (contains the energy for the state). In the second script the vibrational wave functions for the two states and the vibrationally averaged TDMs are now computed using the VIBROT program. We will retain the RASSCF outputs in the scratch directory to check the wave function. It is always dangerous to assume that the wave functions will be correct in a CASSCF calculation. Different problems such as root flippings or incorrect orbitals rotating into the active space are not uncommon. Also, it is always necessary to control that the CASSCF calculation has converged. The first script (Korn shell) is:
#!/bin/ksh
#
# perform some initializations
#
export Project='C2'
export WorkDir=/temp/$LOGNAME/$Project
export Home=/u/$LOGNAME/$Project
echo "No log" > current.log
trap 'cat current.log ; exit 1' ERR
mkdir $WorkDir
cd $WorkDir
#
# Loop over the geometries and generate input for vibrot
#
list="1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 5.0 10.0"
scf='yes'
print "Sigma" > e1.list
print "Pi" > e2.list
for geom in $list
do
# run seward
print "Dist $geom" >> geom.list
sed e "s/#/$geom/" $Home/$Project.seward.input > seward.input
molcas seward.input > current.log
# optionally run scf, motra, guga and mrci to obtain good starting orbitals
if [ "$scf" = 'yes' ]
then
scf='no'
molcas $Home/$Project.scf.input > current.log
molcas $Home/$Project.motra.input > current.log
molcas $Home/$Project.guga.input > current.log
molcas $Home/$Project.mrci.input > current.log
cp $Project.CiOrb $Project.RasOrb1
cp $Project.CiOrb $Project.RasOrb2
fi
# rasscf wavefunction for 1Sg+
ln fs $Project.Job001 JOBIPH
ln fs $Project.RasOrb1 INPORB
molcas $Home/$Project.rasscf1.input > current.log
cat current.log >> rasscf1.log
cat current.log  grep i 'average ci' >> e1.list
cp $Project.RasOrb $Project.RasOrb1
rm f JOBIPH INPORB
# rasscf wavefunction for 1Pu
ln fs $Project.Job002 JOBIPH
ln fs $Project.RasOrb2 INPORB
molcas $Home/$Project.rasscf2.input > current.log
cat current.log >> rasscf2.log
cat current.log  grep i 'average ci' >> e2.list
cp $Project.RasOrb $Project.RasOrb2
rm f JOBIPH INPORB
# rassi to obtain transition
ln fs $Project.Job001 JOB001
ln fs $Project.Job002 JOB002
molcas $Home/$Project.rassi.input > current.log
awk f $Home/tdm.awk current.log >> tdm.list
rm f JOB001 JOB002
#
done
#
# Finished so clean up the files.
#
print "Calculation finished" >&2
cd 
rm $WorkDir/molcas.temp*
#rm r $WorkDir
exit 0
In a second script we will compute the vibrational wave functions
#!/bin/ksh
#
# perform some initializations
#
export Project='C2'
export WorkDir=/temp/$LOGNAME/$Project
export Home=/u/$LOGNAME/$Project
echo "No log" > current.log
trap 'cat current.log ; exit 1' ERR
mkdir $WorkDir
cd $WorkDir
#
# Build vibrot input
#
cp e1.list $Home
cp e2.list $Home
cp geom.list $Home
cp tdm.list $Home
#
cat e1.list geom.list  awk f $Home/wfn.awk > vibrot1.input
cat e2.list geom.list  awk f $Home/wfn.awk > vibrot2.input
cat tdm.list geom.list  awk f $Home/tmc.awk > vibrot3.input
#
ln fs $Project.VibWvs1 VIBWVS
molcas vibrot1.input > current.log
cat current.log
rm f VIBWVS
#
ln fs $Project.VibWvs2 VIBWVS
molcas vibrot2.input > current.log
cat current.log
rm f VIBWVS
#
ln fs $Project.VibWvs1 VIBWVS1
ln fs $Project.VibWvs2 VIBWVS2
molcas vibrot3.input > current.log
cat current.log
rm f VIBWVS1 VIBWVS2
#
# Finished so clean up the files.
#
print "Calculation finished" >&2
cd 
rm $WorkDir/molcas.temp*
#rm r $WorkDir
exit 0
The input for the first part of the calculations include the SEWARD, SCF, MOTRA, GUGA, and MRCI inputs:
&SEWARD &END
Title
C2
Pkthre
1.0D11
Symmetry
X Y Z
Basis set
C.ANOS...3s2p.
C .00000000 .00000000 1.4
End of basis
End of input
&SCF &END
Title
C2
ITERATIONS
40
Occupied
2 1 1 0 2 0 0 0
End of input
&MOTRA &END
Title
C2 molecule
Frozen
1 0 0 0 1 0 0 0
LumOrb
End of input
&GUGA &END
Title
C2 molecule
Electrons
8
Spin
1
Inactive
1 1 1 0 1 0 0 0
Active
0 0 0 0 0 0 0 0
CiAll
1
End of Input
&MRCI &END
Title
C2 molecule
SDCI
End of input
We are going to use a small ANO [3s2p] basis set because our purpose it is not to obtain an extreme accuracy. In the SEWARD input the sign '#' will be replaced by the right distance using the 'sed' command. In the MOTRA input we have frozen the two core orbitals in the molecule, which will be recognized by the MRCI program. The GUGA input defines the reference space of configurations for the subsequent MRCI or ACPF calculation. In this case the valence orbitals are doubly occupied and there is only one reference configuration (they are included as inactive). We thus use one single configuration to perform the SDCI calculation and obtain the initial set of orbitals for the CASSCF calculation.
The lowest state in C_{2} is the result of the electronic configuration [core](2)^{2} (2)^{2} (1)^{4}. Only one electronic state is obtained from this configuration. The configuration (1)^{3} (3)^{1} is close in energy and generates two possibilities, one and one state. The former is the lowest state of the Swan bands, and was thought to be the ground state of the molecule. Transitions to the state are known as the Phillips band and this is the state we are going to compute. We have the possibility to compute the state in symmetry or ( MOLCAS symmetry groups 2 and 3, respectively ) in the D_{2h} group, because both represent the degenerate symmetry in .
The RASSCF input file to compute the two states are:
&RASSCF &END
Title
C2 1Sigmag+ state.
Nactel
4 0 0
Spin
1
Symmetry
1
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
*Average
*2 2 3 6 7
OutOrbitals
Natural
1
Iter
50,25
Lumorb
End of input
&RASSCF &END
Title
C2 1Piu state.
Nactel
4 0 0
Spin
1
Symmetry
2
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
Average
2 2 3 6 7
OutOrbitals
Natural
1
Iter
50,25
Lumorb
End of input
We can skip the SUPSym option because our basis set contains only s,p functions and no undesired rotations can happen. Symmetries and on one hand and and on the other are averaged. Notice that to obtain natural orbitals we have used keyword OUTOrbitals instead of the old RASREAD program. In addition, we need the RASSI input:
&RASSI &END
NrOfJobiphs
2 1 1
1
1
End of input
The VIBROT inputs to compute the vibrationalrotational analysis and spectroscopic constants of the state should be:
&VIBROT &END
RoVibrational spectrum
Title
VibRot spectrum for C2. 1Sigmag+
Atoms
0 C 0 C
Grid
400
Range
2.0 10.0
Vibrations
3
Rotations
0 4
Orbital
0
Potential
2.2 75.42310136
...
End of input
Under the keyword POTEntial the bond distance and potential energy (both in au) of the corresponding state must be included. In this case we are going to compute three vibrational quanta and four rotational quantum numbers. For the state, the keyword ORBItal must be set to one, corresponding to the orbital angular momentum of the computed state. VIBROT fits the potential curve to an analytical curve using splines. The rovibrational Schrödinger equation is then solved numerically (using Numerov's method) for one vibrational state at a time and for the specified number of rotational quantum numbers. File VIBWVS will contain the corresponding wave function for further use.
Just to give some of the results obtained, the spectroscopic constants for the state were:
Re(a) 1.4461
De(ev) 3.1088
D0(ev) 3.0305
we(cm1) .126981E+04
wexe(cm1) .130944E+02
weye(cm1) .105159E+01
Be(cm1) .134383E+01
Alphae(cm1) .172923E01
Gammae(cm1) .102756E02
Dele(cm1) .583528E05
Betae(cm1) .474317E06
and for the state:
Re(a) 1.3683
De(ev) 2.6829
D0(ev) 2.5980
we(cm1) .137586E+04
wexe(cm1) .144287E+02
weye(cm1) .292996E+01
Be(cm1) .149777E+01
Alphae(cm1) .328764E01
Gammae(cm1) .186996E02
Dele(cm1) .687090E05
Betae(cm1) .259311E06
To compute vibrationally averaged TDMs the VIBROT input must be:
&VIBROT &END
Transition moments
Observable
Transition dipole moment
2.2 0.412805
...
End of input
Keyword OBSErvable indicates the start of input for radial functions of observables other than the energy. In the present case the vibrationalrotational matrix elements of the transition dipole moment function will be generated. The values of the bond distance and the TDM at each distance must be then included in the input. VIBROT also requires the VIBWVS1 and VIBWVS2 files containing the vibrational wave functions of the involved electronic states. The results obtained contain matrix elements, transition moments over vibrational wave functions, and the lifetimes of the transition among all the computed vibrationalrotational states. The radiative lifetime of a vibrational level depends on the sum of the transition probabilities to all lower vibrational levels in all lower electronic states. If rotational effects are neglected, the lifetime () can be written as
(10.1) 
where v' and v'' are the vibrational levels of the lower and upper electronic state and A_{v'v''} is the Einstein A coefficient (ns^{1}) computed as
(10.2) 
is the energy difference (au) and TDM_{v'v''} the transition dipole moment (au) of the transition.
For instance, for rotational states zero of the state and one of the state:
Rotational quantum number for state 1: 0, for state 2: 1

~
Overlap matrix for vibrational wave functions for state number 1
1 1 .307535 2 1 .000000 2 2 .425936 3 1 .000000 3 2 .000000 3 3 .485199
~
Overlap matrix for vibrational wave functions for state number 2
1 1 .279631 2 1 .000000 2 2 .377566 3 1 .000000 3 2 .000000 3 3 .429572
~
Overlap matrix for state 1 and state 2 functions
.731192 .617781 .280533
.547717 .304345 .650599
.342048 .502089 .048727
~
Transition moments over vibrational wave functions (atomic units)
.286286 .236123 .085294
.218633 .096088 .240856
.125949 .183429 .005284
~
Energy differences for vibrational wave functions(atomic units)
1 1 .015897 2 1 .010246 2 2 .016427 3 1 .004758 3 2 .010939 3 3 .017108
~
Contributions to inverse lifetimes (ns1)
No degeneracy factor is included in these values.
1 1 .000007 2 1 .000001 2 2 .000001 3 1 .000000 3 2 .000001 3 3 .000000
~
Lifetimes (in nano seconds)
v tau
1 122090.44
2 68160.26
3 56017.08
Probably the most important caution when using the VIBROT program in diatomic molecules is that the number of vibrational states to compute and the accuracy obtained depends strongly on the computed surface. In the present case we compute all the curves to the dissociation limit. In other cases, the program will complain if we try to compute states which lie at energies above those obtained in the calculation of the curve.
This section is a brief comment on a complex situation in a diatomic molecule such as Ni_{2}. Our purpose is to compute the ground state of this molecule. An explanation of how to calculate it accurately can be found in ref. [209]. However we will concentrate on computing the electronic states at the CASSCF level.
The nickel atom has two close lowlying configurations 3d^{8}4s^{2} and 3d^{9}4s^{1}. The combination of two neutral Ni atoms leads to a Ni_{2} dimer whose ground state has been somewhat controversial. For our purposes we commence with the assumption that it is one of the states derived from 3d^{9}4s^{1} Ni atoms, with a single bond between the 4s orbitals, little 3d involvement, and the holes localized in the orbitals. Therefore, we compute the states resulting from two holes on orbitals: states.
We shall not go through the procedure leading to the different electronic states that can arise from these electronic configurations, but refer to the Herzberg book on diatomic molecules [211] for details. In we have three possible configurations with two holes, since the orbitals can be either gerade (g) or ungerade (u): ()^{2}, ()^{1}()^{1}, or ()^{2}. The latter situation corresponds to nonequivalent electrons while the other two to equivalent electrons. Carrying through the analysis we obtain the following electronic states:
()^{2} : , ,
()^{2} : , ,
()^{1}()^{1}: , , ,
, ,
In all there are thus 12 different electronic states.
Next, we need to classify these electronic states in the lower symmetry D_{2h}, in which MOLCAS works. This is done in Table , which relates the symmetry in to that of D_{2h}. Since we have only , , and states here, the D_{2h} symmetries will be only A_{g}, A_{u}, B_{1g}, and B_{1u}. The table above can now be rewritten in D_{2h}:
()^{2} : ( +
),
,
()^{2} : ( +
),
,
()^{1}()^{1}: ( +
), ( +
),
, ,
,
or, if we rearrange the table after the D_{2h} symmetries:
: ()^{2}, ()^{2},
()^{2}, ()^{2}
: ()^{1}()^{1},
()^{1}()^{1}
: ()^{2}, ()^{2}
: ()^{1}()^{1},
()^{1}()^{1}
: ()^{1}()^{1},
()^{1}()^{1}
: ()^{2}, ()^{2}
: ()^{1}()^{1},
()^{1}()^{1}
It is not necessary to compute all the states because some of them (the states) have degenerate components. It is both possible to make single state calculations looking for the lowest energy state of each symmetry or stateaverage calculations in each of the symmetries. The identification of the states can be somewhat difficult. For instance, once we have computed one state it can be a or a state. In this case the simplest solution is to compare the obtained energy to that of the degenerate component in B_{1g} symmetry, which must be equal to the energy of the state computed in A_{g} symmetry. Other situations can be more complicated and require a detailed analysis of the wave function.
It is important to have clean dorbitals and the SUPSym keyword may be needed to separate and (and if gtype functions are used in the basis set) orbitals in symmetry 1 (A_{g}). The AVERage keyword is not needed here because the and orbitals have the same occupation for and states.
Finally, when states of different multiplicities are close in energy, the spinorbit coupling which mix the different states should be included. The CASPT2 study of the Ni_{2} molecule in reference [209], after considering all the mentioned effects determined that the ground state of the molecule is a 0_{g}^{+} state, a mixture of the and electronic states. For a review of the spinorbit coupling and other important coupling effects see reference [214].
There are a large number of symmetry point groups in which MOLCAS cannot directly work. Although unusual in organic chemistry, some of them can be easily found in inorganic compounds. Systems belonging for instance to threefold groups such as C_{3v}, D_{3h}, or D_{6h}, or to groups such O_{h} or D_{4h} must be computed using lower symmetry point groups. The consequence is, as in linear molecules, that orbitals and states belonging to different representations in the actual groups, belong to the same representation in the lower symmetry case, and vice versa. In the RASSCF program it is possible to prevent the orbital and configurational mixing caused by the first situation. The CLEAnup and SUPSymmetry keywords can be used in a careful, and somewhat tedious, way. The right symmetry behaviour of the RASSCF wave function is then assured. It is sometimes not a trivial task to identify the symmetry of the orbitals in the higher symmetry representation and which coefficients must vanish. In many situations the ground state wave function keeps the right symmetry (at least within the printing accuracy) and helps to identify the orbitals and coefficients. It is more frequent that the mixing happens for excited states.
The reverse situation, that is, that orbitals (normally degenerated) which belong to the same symmetry representation in the higher symmetry groups belong to different representations in the lower symmetry groups cannot be solved by the present implementation of the RASSCF program. The AVERage keyword, which performs this task in the linear molecules, is not prepared to do the same in nonlinear systems. Provided that the symmetry problems mentioned in the previous paragraph are treated in the proper way and the trial orbitals have the right symmetry, the RASSCF code behaves properly.
There is a important final precaution concerning the high symmetry systems: the geometry of the molecule must be of the right symmetry. Any deviation will cause severe mixings. Figure contains the SEWARD input for the magnesium porphirin molecule. This is a D_{4h} system which must be computed D_{2h} in MOLCAS.
For instance, the x and y coordinates of atoms C1 and C5 are interchanged with equal values in D_{4h} symmetry. Both atoms must appear in the SEWARD input because they are not independent by symmetry in the D_{2h} symmetry in which MOLCAS is going to work. Any deviation of the values, for instance to put the y coordinate to 0.681879 Å in C1 and the x to 0.681816 Å in C5 and similar deviations for the other coordinates, will lead to severe symmetry mixtures. This must be taken into account when geometry data are obtained from other program outputs or data bases.
The situation can be more complex for some threefold point groups such as D_{3h} or C_{3v}. In these cases it is not possible to input in the exact cartesian geometry, which depends on trigonometric relations and relies on the numerical precision of the coordinates entry. It is necessary then to use in the SEWARD input as much precision as possible and check on the distance matrix of the SEWARD output if the symmetry of the system has been kept at least within the output printing criteria.