%! tutorials.tex $ this file belongs to the Molcas repository $

\chapter{Problem Based Tutorials}

\index{Problem Based Tutorials}

\section{Electronic Energy at Fixed Nuclear Geometry}

The \molcas\ \molcasversion\ suite of Quantum Chemical programs is modular in
design, and a desired calculation is achieved by executing a list of
\molcas\ program modules in succession, occasionally manipulating
the program information files. If the information files from a previous
calculation are saved, then a subsequent calculation need not recompute
them.  This is dependent on the correct information being preserved in
the information files for the subsequent calculations. Each module has keywords 
to specify the functions to be carried out, and many modules rely on the
specification of keywords in previous modules.

In the present examples the calculations will be designed by preparing
a single file in which the input for the different programs is presented
sequentially. The initial problem will be to compute an electronic energy
at a fixed geometry of the nuclei, and this will be performed using different
methods and thus requiring different \molcas\ program modules.

First, the proper \molcas\ environment has to be set up which requires that 
following variables must be properly defined, for instance:

\begin{inputlisting}
export MOLCAS=/home/molcas/molcas.\molcasversion
export Project=CH4
export WorkDir=/home/user/tmp
\end{inputlisting}

If not defined, \molcas\ provides default values for the above environment variables:
\begin{itemize}
\item The {MOLCAS} variable will be set to the latest implemented version of the code.

This variable is set directly in the \molcas\ home directory 

\item Project and WorkDir have the default values None and \$PWD, respectively.

It is very important that the molcas driver, called by command \command{molcas},
and built during the installation of the code, is included in the \$PATH.
\end{itemize}

The first run involves a calculation of the SCF energy of the methane
(CH$_4$) molecule. Three programs should be used: \program{GATEWAY} to specify 
information about the system, \program{SEWARD} to compute
and store the one- and two-electron integrals, and \program{SCF} to obtain
the Hartree-Fock SCF wave function and energy.  

The three \molcas\ programs to 
be used leads to three major entries in the input file: \program{GATEWAY}, \program{SEWARD}, and \program{SCF}.
The \program{GATEWAY} program contains the nuclear geometry in cartesian
coordinates and the label for the one-electron basis set.
The keyword \keyword{coord} allows automatic insertion of \program{GATEWAY} input from a standard
file containing the cartesian coordinates in Angstrom which can be generated by
programs like \program {LUSCUS} or \program{MOLDEN}). 
No symmetry is being considered so the keyword \keyword{group=C1} is used to force the program not 
to look for symmetry in the CH$_4$ molecule, and ,thus, input for \program{SEWARD} is not required. 
In closed-shell cases, like CH$_4$, input for \program{SCF} is not required. All the input
files discussed here can be found at {$\$MOLCAS/doc/samples/problem\_based\_tutorials$}, including the file
\file{SCF.energy.CH4} described below.

%%%To_extract{/doc/samples/problem_based_tutorials/SCF.energy.CH4.input}
\begin{inputlisting}
*SCF energy for CH4 at a fixed nuclear geometry.
*File: SCF.energy.CH4
*
&GATEWAY
 Title = CH4 molecule
 coord  = CH4.xyz 
 basis = STO-3G 
 group = C1

&SEWARD                                                                                                                                                                             
&SCF
\end{inputlisting}
%%%To_extract

where the content of the \file{CH4.xyz} file is: 

%%%To_extract{/doc/samples/problem_based_tutorials/CH4.xyz}
\begin{inputlisting}
5
distorted CH4 coordinates in Angstroms 
C    0.000000     0.000000     0.000000
H    0.000000     0.000000     1.050000
H    1.037090     0.000000    -0.366667
H   -0.542115    -0.938971    -0.383333
H   -0.565685     0.979796    -0.400000
\end{inputlisting}
%%%To_extract

To run \molcas\ , simply execute the command

\begin{inputlisting}
molcas SCF.energy.CH4.input > SCF.energy.CH4.log 2 > SCF.energy.CH4.err
\end{inputlisting}
where the main output is stored in file \file{SCF.energy.CH4.log} 

or

\begin{inputlisting}
molcas -f SCF.energy.CH4.input
\end{inputlisting}
where the main output is stored in \file{SCF.energy.CH4.log}, and the default error file in \file{SCF.energy.CH4.err}. 

The most relevant information is contained in the output file, where the \program{GATEWAY} program 
information describing the nuclear geometry, molecular symmetry, and the data 
regarding the one-electron basis sets and the calculation of one- and 
two-electron integrals, as described in section~\ref{TUT:sec:seward}. Next,
comes the output of program \program{SCF} with information of the electronic
energy, wave function, and the Hartree-Fock (HF) molecular orbitals 
(see section~\ref{TUT:sec:scf}).

Files containing intermediate information, integrals, orbitals, etc, will be 
kept in the {\$WorkDir} directory for further use. For instance, files
\file{\$Project.OneInt} and \file{\$Project.OrdInt} contain the one- and
two-electron integrals stored in binary format. File \file{\$Project.ScfOrb}
stores the HF molecular orbitals in ASCII format, and 
\file{\$Project.RunFile} is a communication file between programs. All these 
files can be used later for more advanced calculations avoiding a 
repeat of certain calculations.

There are graphical utilities that can be used for the analysis of the
results. By default, \molcas\ generates files which can be read with the 
\program{MOLDEN} program and are found in the {\$WorkDir} including the file{CH4.scf.molden}. 
This file contains information about molecular geometry and molecular orbitals, and requires the use if \textit{Density Mode} in \program{MOLDEN}.
However, \molcas\ has its own graphical tool, program \program{LUSCUS}, which is a viewer based on openGL and allows the visualization of 
molecular geometries, orbitals, densities, and density differences. For 
example, a graphical display of the CH$_4$ molecule can be obtained from a standard coordinate file by the following command:

\begin{inputlisting}
luscus CH4.xyz 
\end{inputlisting}

In order to obtain the information for displaying molecular orbitals and densities,
it is necessary to run the \molcas\ program called \program{GRID\_IT}:

%%%To_extract{/doc/samples/problem_based_tutorials/SCF.energy_grid.CH4.input}

\begin{inputlisting}
*SCF energy for CH4 at a fixed nuclear geometry plus a grid for visualization.
*File: SCF.energy_grid.CH4
*
&GATEWAY
 Title = CH4 molecule
 coord = CH4.xyz 
 basis = STO-3G 
 Group = C1

&SEWARD; &SCF                                                                                                                                                                            

&GRID_IT 
 All
\end{inputlisting}
%%%To_extract

Now, execcute the \molcas\ program:

\begin{inputlisting}
molcas SCF.energy_grid.CH4.input -f 
\end{inputlisting}

In the {\$WorkDir} and {\$PWD} directories a new file is generated, \file{CH4.lus} which
contains the information required by the \program{GRID\_IT} input. The file can 
be visualized by \program{LUSCUS} (Open source program, which can be downloaded and 
installed to your Linux, Windows, or MacOS workstation or laptop). By typing the command:

\begin{inputlisting}
luscus CH4.lus
\end{inputlisting}

a window will be opened displaying the molecule and its charge density. By proper
selection of options with the mouse buttons, the shape and size of several molecular orbitals
can be visualized.

\program{GRID\_IT} can also be run separately, if an orbital file is specified in
the input, and the {\$WorkDir} directory is available.

More information can be found in section \ref{UG:sec:gridit}.

As an alternative to running a specific project, the short script provided below can be placed
in the directory {\$MOLCAS/doc/samples/problem\_based\_tutorials} with the name \file{project.sh}.
Simply execute the shell script, \command{project.sh \$Project}, where {\$Project} is the {MOLCAS} input,
and output files, error files, and a {\$WorkDir} directory called {\$Project.work} will be obtained.

%%%To_extract{/doc/samples/problem_based_tutorials/project.sh}
\begin{inputlisting}
#!/bin/bash
                                                                                                                                                                            
export MOLCAS=$PWD
export MOLCAS_DISK=2000
export MOLCAS_MEM=64
export MOLCAS_PRINT=3
                                                                                                                                                                            
export Project=$1
export HomeDir=$MOLCAS/doc/samples/problem_based_tutorials
export WorkDir=$HomeDir/$Project.work
mkdir $WorkDir 2>/dev/null
molcas $HomeDir/$1 >$HomeDir/$Project.log 2>$HomeDir/$Project.err
exit
\end{inputlisting}
%%%To_extract

In order to run a Kohn-Sham density functional calculation, \molcas\ uses the 
same \program{SCF} module, and, therefore, the only change needed are the specification 
of the DFT option and required functional (e.g. B3LYP) in the \program{SCF} input:

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.energy.CH4.input}
\begin{inputlisting}
*DFT energy for CH4 at a fixed nuclear geometry plus a grid for visualization.
*File: DFT.energy.CH4
*
&GATEWAY
 Title = CH4 molecule
 coord = CH4.xyz 
 basis = STO-3G 
 group = C1
&SEWARD
&SCF 
 KSDFT = B3LYP
&GRID_IT 
 All
\end{inputlisting}
%%%To_extract

Similar graphical files can be found in \$WorkDir and \$PWD.

The next step is to obtain  the second-order M{\o}ller--Plesset perturbation (MP2)
energy for methane at the same molecular geometry using the same one-electron
basis set. Program \program{MBPT2} is now used, and it is possible to take 
advantage of having previously computed the proper integrals with \program{SEWARD}
and the reference closed-shell HF wave function with the \program{SCF} program.
In such cases, it is possible to keep the same definitions as before and simply prepare a file 
containing the \program{MBPT2} input and run it using the \command{molcas}
command.

The proper intermediate file will be already in \$WorkDir.
On the other hand, one has to start from scratch, all required inputs should
be placed sequentially in the \file{MP2.energy.CH4} file.
If the decision is to start the project from the beginning,  it is probably a good idea to remove
the entire {\$WorkDir} directory, unless it is known for certain the exact nature of the files contained in this directory.


%%%To_extract{/doc/samples/problem_based_tutorials/MP2.energy.CH4.input}
\begin{inputlisting}
*MP2 energy for CH4 at a fixed nuclear geometry.
*File: MP2.energy.CH4
*
&GATEWAY
 Title = CH4 molecule
 coord = CH4.xyz 
 basis = STO-3G 
 group = C1
&SEWARD 
&SCF
&MBPT2 
 Frozen = 1
\end{inputlisting}
%%%To_extract

In addition to the calculation of a HF wave function, an MP2 calculation has been performed with 
a frozen deepest orbital, the carbon 1s, of CH$_4$. Information about the output
of the \program{MBPT2} program can be found on section~\ref{TUT:sec:mbpt2}.

The \program{SCF} program works by default with closed-shell systems with an
even number of electrons at the Restricted Hartee-Fock (RHF) level. If, 
instead there is a need to use the Unrestricted Hartree Fock (UHF) method, this can be schieved by invoking the
keyword \keyword{UHF}. This is possible for both even and odd electron systems. 
For instance, in a system with an odd number of electrons such as the CH$_3$ radical, with the 
following Cartesian coordinates

%%%To_extract{/doc/samples/problem_based_tutorials/CH3.xyz}
\begin{inputlisting}
4
CH3 coordinates in Angstrom 
C    0.000000     0.000000     0.000000
H    0.000000     0.000000     1.050000
H    1.037090     0.000000    -0.366667
H   -0.542115    -0.938971    -0.383333
\end{inputlisting}
%%%To_extract

the input to run an open-shell UHF calculation is easily obtained:

%%%To_extract{/doc/samples/problem_based_tutorials/SCF.energy_UHF.CH3.input}
\begin{inputlisting}
*SCF/UHF energy for CH3 at a fixed nuclear geometry
*File: SCF.energy_UHF.CH3
*
&GATEWAY
 Title = CH3 molecule
 coord = CH3.xyz 
 basis = STO-3G 
 group = C1
&SEWARD
&SCF 
 UHF
\end{inputlisting}
%%%To_extract

If the system is charged, this must be indicated in the 
\program{SCF} input, for example,  by computing the cation of the CH$_4$ molecule 
at the UHF level:

%%%To_extract{/doc/samples/problem_based_tutorials/SCF.energy_UHF.CH4plus.input}
\begin{inputlisting}
*SCF/UHF energy for CH4+ at a fixed nuclear geometry
*File: SCF.energy_UHF.CH4plus
*
&GATEWAY
 Title = CH4+ molecule
 coord = CH4.xyz 
 basis = STO-3G 
 group = c1
&SEWARD
&SCF
 UHF
 Charge = +1
\end{inputlisting}
%%%To_extract

The Kohn-Sham DFT calculation can be also run using the UHF algorithm:

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.energy.CH4plus.input}
\begin{inputlisting}
*DFT/UHF energy for CH4+ at a fixed nuclear geometry
*File: DFT.energy.CH4plus
*
&GATEWAY
 Title = CH4+ molecule
 coord = CH4.xyz 
 basis = STO-3G 
 group = C1
&SEWARD
&SCF 
 KSDFT = B3LYP
 UHF
 Charge = +1
\end{inputlisting}
%%%To_extract

For the UHF and UHF/DFT methods it is also possible to specify 
$\alpha$ and $\beta$ orbital occupations in two ways. 
\begin{enumerate}
\item First, the keyword \keyword{ZSPIn} can be invoked in the \program{SCF} program, which represents the 
difference between the number of $\alpha$ and $\beta$ electrons. 

For example, setting the keyword to 2 forces the program to converge to a result with two more $\alpha$ than $\beta$ electrons.

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.energy_zspin.CH4.input}
\begin{inputlisting}
*DFT/UHF energy for different electronic occupation in CH4 at a fixed nuclear geometry
*File: DFT.energy_zspin.CH4
*
&GATEWAY
 Title = CH4 molecule 
 coord = CH4.xyz 
 basis = STO-3G 
 group = c1
&SEWARD
&SCF
 Title = CH4 molecule zspin 2
 UHF; ZSPIN =  2
 KSDFT =  B3LYP
\end{inputlisting}
%%%To_extract

The final occupations in the output will show six $\alpha$ and four $\beta$ orbitals.

\item Alternatively, instead of \keyword{ZSPIn}, it is possible to specify 
occupation numbers with keyword \keyword{Occupation} at the beginning of the SCF calculation. 

This requires an additional input line containing the occupied $\alpha$ orbitals (e.g. 6 in this case), and a second line 
with the $\beta$ orbitals (e.g. 4 in this case). Sometimes, SCF convergence may be improved by using this option.
\end{enumerate}

Different sets of methods use other \molcas\ modules. For example, to perform a Complete
Active Space (CAS) SCF calculation, the \program{RASSCF} program has to be used. This
module requires starting trial orbitals, which can be obtained from a previous SCF
calculation or, automatically, from the \program{SEWARD} program which provides trial orbitals by
using a model Fock operator. 

Recommended keywords are 
\begin{itemize}
\item \keyword{Nactel} defines the total number of active 
electrons, holes in Ras1, and particles in Ras3, respectively.  The last two values 
are only for RASSCF-type calculations. 
\item \keyword{Inactive} indicates the number of inactive orbitals where the occupation is always 2 in the CASSCF reference, and 
\item \keyword{Ras2} defines the number of active orbitals. 

By default, the wave function for the lowest state corresponds to the symmetry with spin multiplicity of 1.
Most of the input may not be necessary, if one has prepared and linked an INPORB file with the different orbital types defined by 
a program like \program{LUSCUS}.
\end{itemize}

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.energy.CH4.input}
\begin{inputlisting}
*CASSCF energy for CH4 at a fixed nuclear geometry
*File: CASSCF.energy.CH4
*
&GATEWAY
 coord = CH4.xyz
 basis = STO-3G
 group = C1
&SEWARD
&RASSCF
 Title = CH4 molecule
 Spin = 1; Nactel = 8 0 0; Inactive = 1; Ras2 = 8
&GRID_IT
 All
\end{inputlisting}
%%%To_extract

In this case, the lowest singlet state (i.e. the ground dstate) is computed, since this is a 
closed-shell situation with an active space of eight electrons in eight orbitals and 
with an inactive C 1s orbital, the lowest orbital of the CH$_4$ molecule. This is a CASSCF example in which all the valence 
orbitals and electrons (C 2s, C 2p and 4 x H 1s) are included 
in the active space and allows complete dissociation into
atoms. If this is not the goal, then the three almost degenerate 
highest energy occupied orbitals and the corresponding antibonding unoccupied orbitalsmust be active, leading to 
a 6 in 6 active space.

Using the CASSCF wave function as a reference, it is possible to perform a second-order 
perturbative, CASPT2, correction to the electronic energy by employing the 
\program{CASPT2} program. If all previously calculated files are retained in the 
\$WorkDir directory, in particular, integral files (\file{CH4.OneInt},\file{CH4.OrdInt}), 
the CASSCF wave function information file (\file{CH4.JobIph}), and communication file \file{CH4.RunFile}), it will not be 
necessary to re-run programs \program{SEWARD}, and \program{RASSCF}. In this case
case, it is enough to prepare a file containing input only for the \program{CASPT2} program followed be execution.
Here, however, for the sake of completness, input to all \molcas\ moddules is provided:

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.energy.CH4.input}
\begin{inputlisting}
*CASPT2 energy for CH4 at a fixed nuclear geometry
*File: CASPT2.energy.CH4
*
&GATEWAY
 coord = CH4.xyz; basis = STO-3G; group = C1
&SEWARD
&RASSCF
   LumOrb
   Spin = 1; Nactel = 8 0 0; Inactive = 1; Ras2 = 8
&CASPT2
 Multistate = 1 1
\end{inputlisting}
%%%To_extract

In most of casesi, the Hartree-Fock orbitals will be a better choice as starting orbitals. 
In that case, the \program{RASSCF} input has to include keyword \keyword{LumOrb} to read 
from any external source of orbitals other than those generated by the \program{SEWARD} program. 
By modifying input to the \program{SCF} program, it is possible to generate 
alternative trial orbitals for the \program{RASSCF} program. Since a new set of trial orbitals is used,
the input to the \program{RASSCF} program is also changed. Now, the number of
active orbitals, as well as the number of active electrons, are 6. 

The two lowest orbitals (\keyword{Inactive} 2) are excluded from the active space
and one other orbital is placed in the secondary space.
If the previous (8,8) full valence space was used,
the \program{CASPT2} program would not be able to include more electronic correlation energy,
considering that the calculation involves a minimal basis set.
The input for the \program{CASPT2} program includes a frozen C 1s orbital, the lowest orbital
in the CH$_4$ molecule.

The charge and multiplicity of our wave function can be changed by computing the
CH$_4^+$ cation with the same methods. The \program{RASSCF} program defines
the character of the problem by specifying the number of electrons, the spin multiplicity, and the spatial
symmetry. In the example below, there is one less electron giving rise to doublet multiplicity:

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.energy.CH4plus.input}
\begin{inputlisting}
*CASSCF energy for CH4+ at a fixed nuclear geometry
*File: CASSCF.energy.CH4plus
*
&GATEWAY
 Title = CH4+ molecule 
 coord = CH4.xyz; basis = STO-3G; Group = C1
&SEWARD
&RASSCF 
   Spin = 2; Nactel = 7 0 0; Inactive = 1; Ras2 = 8
\end{inputlisting}
%%%To_extract

No further modification is needed to the \program{CASPT2} input:

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.energy.CH4plus.input}
\begin{inputlisting}
*CASPT2 energy for CH4+ at a fixed nuclear geometry
*File: CASPT2.energy.CH4plus
*
&GATEWAY
 coord = CH4.xyz; basis = STO-3G; group = C1
&SEWARD
&RASSCF
   Title = CH4+ molecule
   Spin = 2; Nactel = 1 0 0; Inactive = 4; Ras2 = 1
&CASPT2
\end{inputlisting}
%%%To_extract

A somewhat more sophisticated calculation can be performed at the
Restricted Active Space (RAS) SCF level. In such a situation, the level of excitation
in the CI expansion can be controlled by restricting the number of holes
and particles present in certain orbitals. 

%%%To_extract{/doc/samples/problem_based_tutorials/RASSCF.energy.CH4.input}
\begin{inputlisting}
*RASSCF energy for CH4 at a fixed nuclear geometry
*File: RASSCF.energy.CH4
*
&GATEWAY
 coord = CH4.xyz; basis = STO-3G; group = C1
&SEWARD
&RASSCF
   Title = CH4 molecule
   Spin = 1; Nactel = 8 1 1
   Inactive = 1; Ras1 = 1; Ras2 = 6; Ras3 = 1
\end{inputlisting}
%%%To_extract

In particular, the previous calculation includes one orbital within the Ras1
space and one orbital within the Ras3 space. One hole (single excitation) at
maximum is allowed from Ras1 to Ras2 or Ras3, while a maximum of one particle
is allowed in Ras3, derived from either Ras1 or Ras2. Within Ras2, all types
of orbital occupations are allowed. The RASSCF wave functions can be used 
as reference for multiconfigurational perturbation theory (RASPT2), but
this approach has not been as extensively tested as CASPT2, and, so experience is
still somewhat limited.

\molcas\ also has the possibility of computing electronic energies at 
different CI levels by using the \program{MRCI} program. The input provided below involves
a Singles and Doubles Configuration Interaction  (SDCI) calculation on the CH$_4$ molecule.
To set up the calculations, program \program{MOTRA} which transforms
the integrals to molecular basis, and program \program{GUGA} which computes the
coupling coefficients, must be run before the \program{MRCI} program.
In \program{MOTRA} the reference orbitals are specifiedi, and those employed 
here are from an HF \program{SCF} calculation including frozen orbitals. In \program{GUGA}
the reference for the CI calculation is described by the number of correlated electrons,
the spatial and spin symmetry, the inactive orbitals always occupation 2 in
the reference space, and the type of CI expansion. 

%%%To_extract{/doc/samples/problem_based_tutorials/SDCI.energy.CH4.input}
\begin{inputlisting}
*SDCI energy for CH4 at a fixed nuclear geometry
*File: SDCI.energy.CH4
*
&GATEWAY
 coord = CH4.xyz; basis = STO-3G; group = c1
&SEWARD
&SCF
&MOTRA
 Lumorb
 Frozen= 1
&GUGA 
 Electrons = 8
 Spin = 1
 Inactive= 4
 Active= 0
 Ciall= 1
&MRCI 
 SDCI
\end{inputlisting}
%%%To_extract

To use reference orbitals from a previous CASSCF calculation, the
\program{RASSCF} program will have to be run before the \program{MOTRA}
module. Also, if the spatial or spin symmetry are changed for the CI
calculation, the modifications will be introduced in the input to \program{GUGA} program.
Many alternatives are possible for performing an MRCI calculation as shown in the next example below,
in which the reference space to perform the CI is multiconfigurational:

%%%To_extract{/doc/samples/problem_based_tutorials/MRCI.energy.CH4.input}
\begin{inputlisting}
*MRCI energy for CH4 at a fixed nuclear geometry
*File: MRCI.energy.CH4
*
&GATEWAY
 Title = CH4 molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
&SEWARD
&SCF
&RASSCF
 LumOrb
 Spin= 1; Nactel= 6 0 0; Inactive= 2; Ras2= 6
&MOTRA 
 Lumorb
 Frozen= 1
&GUGA
 Electrons= 8
 Spin= 1
 Inactive= 2
 Active= 3
 Ciall= 1
&MRCI
 SDCI
\end{inputlisting}
%%%To_extract

The \program{MRCI} program also allows the calculation of electronic energies using the
ACPF method. Another \molcas\ program, \program{CPF}, offers the possibility to 
use the CPF, MCPF, and ACPF methods with a single reference function. The 
required input is quite similar to that for the \program{MRCI} program:

%%%To_extract{/doc/samples/problem_based_tutorials/CPF.energy.CH4.input}
\begin{inputlisting}
*CPF energy for CH4 at a fixed nuclear geometry
*File: CPF.energy.CH4
*
&GATEWAY
 Title= CH4 molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
&SEWARD
&SCF 
&MOTRA 
 Lumorb
 Frozen= 1
&GUGA
 Electrons= 8
 Spin = 1
 Inactive = 4
 Active = 0
 Ciall= 1
&CPF
 CPF
End Of Input
\end{inputlisting}
%%%To_extract

Finally, \molcas\ can also perform closed- and open-shell coupled cluster
calculations at the CCSD and CCSD(T) levels. These calculations are controlled by
the \program{CCSDT} program, whose main requirement is that the reference 
function has to be generated with the \program{RASSCF} program. The following input is 
required to obtain a CCSD(T) energy for the CH$_4$ molecule:

%%%To_extract{/doc/samples/problem_based_tutorials/CCSDT.energy.CH4.input}
\begin{inputlisting}
*CCSDT energy for CH4 at a fixed nuclear geometry
*File: CCSDT.energy.CH4
*
&GATEWAY
 Title= CH4 molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
&SEWARD
&RASSCF
 Spin= 1; Nactel= 0 0 0; Inactive= 5; Ras2= 0
 OutOrbitals
 Canonical
&MOTRA
 JobIph
 Frozen= 1
&CCSDT
 CCT
\end{inputlisting}
%%%To_extract

Since this is a closed-shell calculation, the \program{RASSCF} input 
computes a simple RHF wave function with zero active electrons and orbitals using 
keywords \keyword{OutOrbitals} and \keyword{Canonical}. The \program{MOTRA} must
include the keyword \keyword{JobIph} to extract the wave function information
from file \file{JOBIPH} which is automatically generated by \program{RASSCF}. Finally,
the keywork \keyword{CCT} in program \program{CCSDT} leads to the calculation of the
CCSD(T) energy using the default algorithms.

The \program{CCSDT} program in \molcas\ is especially suited to compute open-shell
cases. The input required to obtain the electronic energy of the CH$_4^+$ cation
with the CCSD(T) method is:


%%%To_extract{/doc/samples/problem_based_tutorials/CCSDT.energy.CH4plus.input}
\begin{inputlisting}
*CCSDT energy for CH4+ at a fixed nuclear geometry
*File: CCSDT.energy.CH4plus
*
&GATEWAY
 Title= CH4+ molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
&SEWARD
&RASSCF
 Spin= 2; Nactel= 1 0 0; Inactive= 4; Ras2= 1
 OutOrbitals
 Canonical
&MOTRA
 JobIph
 Frozen= 1
&CCSDT
 CCT
\end{inputlisting}
%%%To_extract

where the \program{RASSCF} program generated the proper Restricted Open-Shell 
Hartree-Fock (ROHF) reference. Different levels of spin adaptation are also available.

If solvent effects are desired, \molcas\ includes two
models: Kirkwood and PCM. Adding a solvent effect to a ground state at HF, DFT, or CASSCF levels,
simply requires the inclusion of the keyword \keyword{RF-input} within the input for the \program{SEWARD} 
which calculates a self-consistend reaction field.

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.energy_solvent.CH4.input}
\begin{inputlisting}
*DFT energy for CH4 in water at a fixed nuclear geometry
*File: DFT.energy_solvent.CH4
*
&GATEWAY
 Title= CH4 molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
 RF-input
   PCM-model; solvent= water
 End of RF-input
&SEWARD
&SCF
KSDFT= B3LYP
\end{inputlisting}
%%%To_extract

Other programs such as \program{CASPT2}, \program{RASSI}, and \program{MOTRA} require that
the reaction field is included as a perturbation with keyword \keyword{RFPErturbation}.
In the next example the correction is added at both the CASSCF and CASPT2 levels.

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.energy_solvent.CH4.input}
\begin{inputlisting}
*CASPT2 energy for CH4 in acetone at a fixed nuclear geometry
*File: CASPT2.energy_solvent.CH4
*
&GATEWAY
 Title= CH4 molecule
 coord = CH4.xyz; basis = STO-3G; group = c1
  RF-input
   PCM-model; solvent= acetone; AAre= 0.2
  End of RF-input
&SEWARD
&RASSCF
  Spin= 1; Nactel= 6 0 0; Inactive= 2; Ras2= 6
&CASPT2
 Frozen= 1
 Multistate= 1 1
 RFPert
\end{inputlisting}
%%%To_extract

%)
Notice that the tesserae of the average area in the PCM model (keyword
has been changed to the value required for acetone by the keyword \keyword{Aare},
while the default is 0.4 \AA$^2$ for water 
\ifmanual
(see section~\ref{UG:sec:rfield}).
More detailed examples can be found in section~\ref{TUT:sec:cavity}. 
\fi

\section{Optimizing geometries}
%: minima, transition states, crossings, and minimum energy paths

It is now useful to explore potential energy surfaces (PES) and optimize the molecular geometry for
specific points along the PES. Different cases are discussed including a way to obtain the optimal geometry
in a minimum energy search, to obtain a transition-state structure connecting different regions of 
the PES, to find the crossing between two PES where the energy becomes degenerate, or to map
the minimum steepest-descent energy path (MEP) from an initial point to the final 
a minimum energy geometry as the PES progresses in a downward manner. 

All these types of searches can be performed either by fully optimizing all 
degrees of freedom of the system or by introducing certain restrictions. \molcas\ \molcasversion\ can perform
geometry optimizations at the SCF (RHF and UHF), DFT (RHF and UHF based), CASSCF (CASSCF and RASSCF) levels of theory, 
where efficient analytical gradients are available and at the CASPT2 and other correlated levels where numerical
gradients are used.

Geometry optimizations require many cycles, in which the electronic energy is estimated at a specific
level of calculation followed by calculation of the gradient of the energy with respect to the geometric
degrees of freedom (DOF). With this information at hand, the program must decide if the molecule is 
already at the final required geometry (i.e. gradient $\sim$ 0 for all
DOF) indicating a minimum in the PES or if the geometry must be modified 
and continue the cycle. The input file should,
therefore, be built in a way that allows a loop over the different programs. 

The general input commands \command{Do while} and \command{Enddo} control the loop 
and program input is inserted within these commands. Instructions for the number of maximum iterations allowed and the type of output required can also be added.
\ifmanual
(see section~\ref{UG:sec:sysvar})
\fi
%The commands \command{Set output file}, which prints output for each iterations and
%in the \$WorkDir directory with the file name Structure.\$iteration.output, and 
%\command{Set maxiter 100}, which sets maximum iterations to one hundred.

All examples previously discussed, use \keyword{COORD} keyword, but it also possible
to use \textit{native format}, where symmetry unique atoms are specified (\keyword{SYMMETRY}) 
and provide generators to construct all atoms in the molecule.

The selected example describes geometry optimization of the water molecule at the SCF RHF level
of calculation:

%%%To_extract{/doc/samples/problem_based_tutorials/Water_distorted.xyz}
\begin{inputlisting}
3
 coordinates for water molecule NOT in equilibrium 
O 0.000000  0.000000  0.000000 
H 0.758602  0.000000  0.504284 
H 0.758602  0.000000 -0.504284 
\end{inputlisting}
%%%To_extract


%%%To_extract{/doc/samples/problem_based_tutorials/SCF.minimum_optimization.H2O.input}
\begin{inputlisting}
*SCF minimum energy optimization for H2O
*File: SCF.minimum_optimization.H2O
*
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1

>>> Do while
 &SEWARD ;&SCF; &SLAPAF
>>> EndDo
\end{inputlisting}
%%%To_extract

The sequence of programs employed includes \program{GATEWAY} which is external to the loop, followed by
\program{SEWARD}, \program{SCF}, and \program{SLAPAF}. \program{SEWARD} 
computes the integrals,  \program{SCF} program computes the RHF energy and wave 
function, and  \program{SLAPAF} will control the calculation of gradients and
estimate if the calculation has already finished or needs to proceed to a new
nuclear geometry for the next iteration. Automatically, a file named 
\file{\$Project.geo.molden} will be generated in \$WorkDir containing all the 
geometric steps contained in the optimization process. \program{MOLDEN} or \program{LUSCUS} can 
then read this file to display the individual molecular geometries which form the optimization cycle.
 
Using another reference wave function can be simply performed by changing the sequence of 
programs. For instance, we can perform an UHF calculation of the H$_2$O$^+$ 
cation:
  
%%%To_extract{/doc/samples/problem_based_tutorials/UHF.minimum_optimization.H2Oplus.input}
\begin{inputlisting}
*UHF minimum energy optimization for H2O+
*File: UHF.minimum_optimization.H2Oplus
*
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1
>> Do while

 &SEWARD 
 &SCF; Title="H2O minimum optimization"; UHF; Charge=1
 &SLAPAF 

>> EndDo
\end{inputlisting}
%%%To_extract

The same procedure can be followed if we pretend to perform a DFT geometry optimization:

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.minimum_optimization.H2O.input}
\begin{inputlisting}
*DFT minimum energy optimization for H2O
*File: DFT.minimum_optimization.H2O
*
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1

>>> Export MOLCAS_MAXITER=100
>>> Do while

 &SEWARD
 &SCF ; Title="H2O minimum optimization"; KSDFT=B3LYP
 &SLAPAF &END

>>> EndDo
\end{inputlisting}
%%%To_extract

Once an energy minimum is found based on the calculation of gradients, it is necessary to
ensure that the geometry really is a minimum energy point. This can be only 
accomplished by computing second derivatives of the energy (i.e. the Hessian). 
\molcas\ can compute analytical Hessians for SCF and single state
CASSCF wave functions. For other methods,  numerical procedures can be used 
to compute the Hessian. Once the Hessian is computed, vibrational
frequencies are calculated, and Statistical Mechanics is used to obtain thermodynamic 
properties. At a true energy minimum, there will be 3N-6 real frequencies 
Program \program{MCKINLEY} computes second derivatives 
of a predefined (SCF or CASSCF) wave function, while \program{MCLR} performs 
the vibrational and statistical analyses. \molcas\ simply requires input for 
the \program{MCKINLEY} program to perform the entire calculation by using keywords
\keyword{Perturbation} and \keyword{Hessian}, while program \program{MCLR} will be 
called automatically but requires no input. 
The full set of calculationsi is included below first a geometry optimization followed by the 
calculation of a Hessian.

%%%To_extract{/doc/samples/problem_based_tutorials/SCF.minimization_plus_Hessian.H2O.input}
\begin{inputlisting}
*SCF minimum energy optimization plus hessian of the water molecule
*File: SCF.minimization_plus_hessian.H2O
*
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1

>>> Export MOLCAS_MAXITER=100
>>> Do while

 &SEWARD
 &SCF; Title="H2O minimum optimization"
 &SLAPAF

>>> EndDo

&MCKINLEY
\end{inputlisting}
%%%To_extract

Note that \program{MCKINLEY} input above is placed after \command{EndDo}, and, therefore,
is external to the looping scheme. Once the geometry optimization at the desired level of theory has finished, the 
Hessian will be computed at the final geometry.
In general, any calculation performed using a \$WorkDir directory where a 
previous geometry optimization has taken place will use the last geomtry calculated 
from that optimization as the input geometry even if \program{SEWARD} input is 
present. To avoid that, the only solution is to remove the communication file 
\file{RUNFILE} where the geometry is stored. Note also, that the frequencies are 
computed in a cartesian basis, and that three translational and three rotational 
frequencies which should be very close to zero are included in the output file.
This is not the case when numerical gradients and Hessians are used. 
In particular, for water at its minimum energy structure three (3N-6) 
real vibrational frequencies. By default, in \$WorkDir a file \file{\$Project.freq.molden}
is generated containing the vibrational frequencies and modes, which can be visualized by \program{MOLDEN}.

A new level of theory, CASSCF, is introduced here which is especially suited for 
geometry optimizations of excited states discussed in the next chapter.
A geometry optimization is performed to illustrate a broader range of possibilities including 
the imposition of a geometric restrain that the HOH angle in water should be constrained to 120$^o$
during the optimization.
This means that only the O-H bond distances be optimized in this partial minimization. 
The restriction is indicated
in in \program{GATEWAY}
by invoking the keyword \keyword{Constraints} and ending with the keyword \keyword{End of Constraints}. 
The names of variables corresponding to geometrical variables in either internal or Cartesian coordinates
that are to be constrained are placed between these two keywords. 
\ifmanual
(see nomenclature in 
section~\ref{UG:sec:definition_of_internal_coordinates})
\fi
In the case of H$_2$O, the H1-O-H2 angle is fixed at 120$^o$, so a variable,
$a$, is first defined with the keywork \keyword(Angle), which relates it to the H1-O1-H2 angle, followed by the second keyword, \keyword{Value}, 
where the variable $a$ is specified as 120$^o$. 
It is not required that the initial geometry is 120$^o$, only that the final result for the calculation
will become 120$^o$.

Note that the \program{RASSCF} program requires initial trial orbitals, and those  
which are automatically generated by \program{SEWARD} are used.  The resulting CASSCF 
wave function includes all valence orbitals and electrons.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.minimum_optimization_restricted.H2O.input}
\begin{inputlisting}
*CASSCF minimum energy optimization of the water molecule with geometrical restrictions
*File: CASSCF.minimum_optimization_restricted.H2O
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1
Constraint
   a = Angle H2 O1 H3
  Value
   a = 90. degree
End of Constraints

>>> Do while
   
 &SEWARD 
 &RASSCF; nActEl=8 0 0; Inactive=1; Ras2=6
 &SLAPAF 

>>> EndDo
\end{inputlisting}
%%%To_extract

Other more flexible ways to impose geometric restrictions involve the specification of which internal
coordinates should remain fixed and which should change. In the next example, 
the bond lengths are forced to remain fixed at their initial distance (here 0.91 \AA), while the
bond angle, having an initial of 81$^\circ$, is optimized.

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.minimum_optimization_restricted.H2O.input}
\begin{inputlisting}
*DFT minimum energy optimization of the angle in the water molecule at fixed bond lengths
*File: DFT.minimum_optimization_restricted.H2O
*
&GATEWAY 
 Title= H2O minimum optimization
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1

>>> EXPORT MOLCAS_MAXITER=100
>>> Do while
                                                                                                                                                                            
 &SEWARD; &SCF; Title="H2O restricted minimum"; KSDFT=B3LYP
 &SLAPAF 
  Internal Coordinates
     b1 = Bond O1 H2
     b2 = Bond O1 H3
     a1 = Angle H2 O1 H3
  Vary
     a1
  Fix
     b1
     b2
  End of Internal

>>> EndDo
\end{inputlisting}
%%%To_extract

In the final output, the two O-H bond lengths remain at the initia values, while the H1-O1=H2 angle is optimized
to a final angle of 112$^o$.

The next step entails the computation of a transition state, a structure connecting different regions of
the potential energy hypersurface, and is a maximum for only one degree of
freedom. The most common saddle points have order one, that is, they are maxima for one of
one displacement and minima for the others. The simplest way to search for a 
transition state in \molcas\ is to add the keyword \keyword{TS} to the 
\program{SLAPAF} input. Keyword \keyword{PRFC} is suggested in order to verify 
the nature of the transition structure. Searching for transition states is, 
however, not an easy task. An illustration of the input required for transition state optimization for water at the DFT level 
is given below:

%%%To_extract{/doc/samples/problem_based_tutorials/Water_TS.xyz}
\begin{inputlisting}
3
water in Transition state in bohr
O1             0.750000        0.000000        0.000000 
H2             1.350000        0.000000        1.550000 
H3             1.350000        0.000000       -1.550000 
\end{inputlisting}
%%%To_extract


%%%To_extract{/doc/samples/problem_based_tutorials/DFT.transition_state.H2O.input}
\begin{inputlisting}
*DFT transition state optimization of the water molecule 
*File: DFT.transition_state.H2O
*
&Gateway
 Coord=Water_TS.xyz
 Basis=ANO-S-VDZ
 Group=C1
>>> Do while

 &SEWARD
 &SCF; Title="H2O TS optimization"; KSDFT=B3LYP
 &SLAPAF ; ITER=20 ; TS

>>> EndDo
\end{inputlisting}
%%%To_extract

The initial coordinates were chosen in units of Bohr, to illustrare that this is the
default case.  The optimal geometry for ground state of water is a structure with C$_{2v}$ symmetry. 
A transition state has been found with a linear HOH angle of 180$^o$. 
In many cases, there may be a clue along the energy pathway for a chemical reaction about the nature of the transition state structure, 
which typically represents an intermediate conformation between reactants and products. 
If this turns out to be the case, it is possible to help the optimization process 
proceed toward an informed guess, by invoking the keyword \keyword{FindTS} in \program{SLAPAF}.
\keyword{FindTS} must to be accompanied with a definition of constrained geometric definitions.
\program{SLAPAF} will guide the optimization of the transition state towards a region in
which the restriction is fulfilled. Once there, the restriction will be released
and a free search of the transition state will be performed. This technique is 
frequently quite effective and makes it possible to find difficult transition 
states or reduce the number of required iterations. Here, an example is provided, in 
which the initial geometry of water is clearly bent, and a trial restraint is imposed
such that the angle for the transition state should be near 180$^o$. The 
final transition state will, however, be obtained without any type of geometrical restriction.

%%%To_extract{/doc/samples/problem_based_tutorials/DFT.transition_state_restricted.H2O.input}
\begin{inputlisting}
*DFT transition state optimization of the water molecule with geometrical restrictions
*File: DFT.transition_state_restricted.H2O
*
&Gateway
 Coord=Water_TS.xyz
 Basis=ANO-S-VDZ
 Group=C1
 Constraints
   a = Angle H2 O1 H3
 Value
   a = 180.0 degree
 End of Constraints

>>> Do while

 &SEWARD
 &SCF; Title="H2O TS optimization"; KSDFT=B3LYP
 &SLAPAF ;FindTS

>>> EndDo
\end{inputlisting}
%%%To_extract

The \program{CASPT2} geometry optimizations are somewhat different because \program{ALASKA}
is not suited to compute \program{CASPT2} analytical gradients. Therefore the \program{ALASKA}
program is automatically substituted by program \program{NUMERICAL\_GRADIENT}, which will take care
of performing numerical gradients. From the user point of view the only requirement is to place
the \program{CASPT2} input after the \program{RASSCF} input.
The CASSCF wave function has of course to be generated in each step before 
performing CASPT2. To compute a numerical gradient can be quite time consuming, 
although it is a task that can be nicely parallelized. In a double-sided 
gradient algorithm like here a total of 6N-12+1 CASPT2 calculations are performed 
each pass of the optimization, where N is the number of atoms.

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.minimum_optimization.H2O.input}
\begin{inputlisting}
*CASPT2 minimum energy optimization for water
*File: CASPT2.minimum_optimization.H2O
*
&GATEWAY 
 coord=Water_distorted.xyz
 basis=ANO-S-MB
 group=C1

>>> Do while

 &SEWARD
 &RASSCF; Title="H2O restricted minimum"; nActEl=8 0 0; Inactive=1; Ras2=6
 &CASPT2; Frozen=1 
 &SLAPAF 

>>> EndDo
\end{inputlisting}
%%%To_extract

The use of spatial symmetry makes the calculations more efficient, although
they may again complicate the preparation of input files. We can repeat the previous \program{CASPT2}
optimization by restricting the molecule to work in the C$_{2v}$ point group, which, by the way,
is the proper symmetry for water in the ground state. The \program{GATEWAY} program (as no symmetry
has been specified) will identify and work with the highest available point group,
C$_{2v}$. Here the molecule is placed with YZ as the molecular plane. By adding
keyword \keyword{Symmetry} containing as elements of symmetry the YZ (symbol X) and YX (symbol Z),
the point group is totally defined and the molecule properly generated. From that point the
calculations will be restricted to use symmetry restrictions. For instance, the molecular
orbitals will be classified in the four elements of symmetry of the group, a$_1$, b$_1$, b$_2$,
and a$_2$, and most of the programs will require to define the selection of the orbitals in
the proper order. The order of the symmetry labels is determined by \program{SEWARD} and must
be checked before proceeding, because from that point the elements of symmetry will be known
by their order in \program{SEWARD}: a$_1$, b$_1$, b$_2$, and a$_2$, for instance, will be
symmetries 1, 2, 3, and 4, respectively. \program{SCF} does not require to specify the
class of orbitals and it can be used as a learning tool.


%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.minimum_optimization_C2v.H2O.input}
\begin{inputlisting}
*CASPT2 minimum energy optimization for water in C2v
*File: CASPT2.minimum_optimization_C2v.H2O
*
 &GATEWAY
Title= H2O caspt2 minimum optimization
Symmetry= X Z
Basis set
O.ANO-S...2s1p.
O        0.000000  0.000000  0.000000 Angstrom
End of basis
Basis set
H.ANO-S...1s.
H1       0.000000  0.758602  0.504284 Angstrom
End of basis

>>> EXPORT MOLCAS_MAXITER=100
>>> Do while

 &SEWARD
 &RASSCF; nActEl=8 0 0; Inactive=1 0 0 0; Ras2=3 1 2 0
 &CASPT2; Frozen=1 0 0 0
 &SLAPAF &END

>>> EndDo
\end{inputlisting}
%%%To_extract

Thanks to symmetry restrictions the number of iterations within \program{NUMERICAL\_GRADIENT}
has been reduced to five instead of seven, because many of the deformations 
are redundant within the C$_{2v}$ symmetry. Also, symmetry considerations are 
important when defining geometrical restrictions 
\ifmanual
(see sections~\ref{UG:sec:definition_of_internal_coordinates}
and \ref{TUT:sec:optim}).
\else
(see online manual).
\fi

\section{Computing excited states}

The calculation of electronic excited states is typically a multiconfigurational problem, and
therefore it should preferably be treated with multiconfigurational methods such as CASSCF and
CASPT2. We can start this section by computing the low-lying electronic states of the
acrolein molecule at the CASSCF level and using a minimal 
basis set. The  standard file with cartesian coordinates is:

%%%To_extract{/doc/samples/problem_based_tutorials/acrolein.xyz}
\begin{inputlisting}
 8
Angstrom
 O      -1.808864   -0.137998    0.000000
 C       1.769114    0.136549    0.000000
 C       0.588145   -0.434423    0.000000
 C      -0.695203    0.361447    0.000000
 H      -0.548852    1.455362    0.000000
 H       0.477859   -1.512556    0.000000
 H       2.688665   -0.434186    0.000000
 H       1.880903    1.213924    0.000000
\end{inputlisting}
%%%To_extract

We shall carry out State-Averaged (SA) CASSCF calculations, in which one single 
set of molecular orbitals is used to compute all the states of a given spatial 
and spin symmetry. The obtained density matrix is the average for all states 
included, although each state will have its own set of optimized CI 
coefficients. Different weights can be considered for each of the states, 
but this should not be used except in very special cases by experts. It is 
better to let the CASPT2 method to handle that. The use of a SA-CASSCF 
procedure has an great advantage. For example, all states in a SA-CASSCF 
calculation are orthogonal to each other, which is not necessarily true for
state specific calculations. Here, we shall include five states of singlet 
character the calculation. As no symmetry is invoked all the states belong by 
default to the first symmetry, including the ground state.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.excited.acrolein.input}
\begin{inputlisting}
*CASSCF SA calculation on five singlet excited states in acrolein
*File: CASSCF.excited.acrolein
*
&GATEWAY
  Title= Acrolein molecule
  coord = acrolein.xyz; basis = STO-3G; group = c1
&SEWARD; &SCF
&RASSCF
  LumOrb
  Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
  CiRoot= 5 5 1
&GRID_IT
  All
\end{inputlisting}
%%%To_extract


We have used as active all the $\pi$ and $\pi^*$ orbitals, two bonding and
two antibonding $\pi$ orbitals with four electrons and in addition the oxygen 
lone pair ($n$). Keyword \keyword{CiRoot} informs the program that we want to 
compute a total of five states, the ground state and the lowest four excited 
states at the CASSCF level and that all of them should have the same weight in 
the average procedure. Once analyzed we find that the calculation has provided,
in this order, the ground state, two $n\to\pi^*$ states, and two $\pi\to\pi^*$ states.
It is convenient to add the \program{GRID\_IT} input in order to be able to use
the \program{LUSCUS} interface for the analysis of the orbitals and the occupations
in the different electronic states. Such an analysis should always be made in 
order to understand the nature of the different excited states.
In order to get a more detailed analysis of the nature of the obtained states it is
also possible to obtain in a graphical way the charge density differences between
to states, typically the difference between the ground and an excited state. The
following example creates five different density files:

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.excited_grid.acrolein.input}
\begin{inputlisting}
*CASSCF SA calculation on five singlet excited states in acrolein
*File: CASSCF.excited_grid.acrolein
*
&GATEWAY
  Title= Acrolein molecule
  coord= acrolein.xyz; basis= STO-3G; group= c1
&SEWARD; &SCF
&RASSCF 
 LumOrb
 Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
 CiRoot= 5 5 1
 OutOrbital
 Natural= 5
&GRID_IT
 FILEORB = $Project.RasOrb.1
 NAME = 1; All
&GRID_IT
 FILEORB = $Project.RasOrb.2
 NAME = 2; All
&GRID_IT
 FILEORB = $Project.RasOrb.3
 NAME = 3; All
&GRID_IT
 FILEORB = $Project.RasOrb.4
 NAME = 4; All
&GRID_IT
 FILEORB = $Project.RasOrb.5
 NAME = 5; All
\end{inputlisting}
%%%To_extract

In \program{GRID\_IT} input we have included all orbitals. It is, however,
possible and in general recommended to restrict the calculation to certain
sets of orbitals. How to do this is described in the input manual for
\program{GRID\_IT}. 

Simple math operations can be performed with grids of the same size, 
for example, \program{LUSCUS} can be used to display the difference 
between two densities. 

CASSCF wave functions are typically good enough, but this is not the case for
electronic energies, and the dynamic correlation effects have to be included,
in particular here with the CASPT2 method. The proper input is prepared, again
including \program{SEWARD} and \program{RASSCF} (unnecessary if they were
computed previously), adding a \program{CASPT2} input with the keyword
\keyword{MultiState} set to 5 1 2 3 4 5. The \program{CASPT2} will perform four
consecutive single-state (SS) CASPT2 calculations using the SA-CASSCF roots computed
by the \program{RASSCF} module. At the end, a multi-state CASPT2 calculation
will be added in which the five SS-CASPT2 roots will be allowed to interact.
The final MS-CASPT2 solutions, unlike the previous SS-CASPT2 states, will be
orthogonal. The \keyword{FROZen} keyword is put here as a reminder. By
default the program leaves the core orbitals frozen.


%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.excited.acrolein.input}
\begin{inputlisting}
*CASPT2 calculation on five singlet excited states in acrolein
*File: CASPT2.excited.acrolein
*
&GATEWAY
 Title= Acrolein molecule
 coord = acrolein.xyz; basis = STO-3G; group= c1
&SEWARD; &SCF
&RASSCF
 Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
 CiRoot= 5 5 1
&GRID_IT
 All
&CASPT2
 Multistate= 5 1 2 3 4 5
 Frozen= 4
\end{inputlisting}
%%%To_extract


Apart from energies and state properties it is quite often necessary to compute
state interaction properties such as transition dipole moments, Einstein coefficients,
and many other. This can be achieved with the \program{RASSI} module, a powerful
program which can be used for many purposes 
\ifmanual
(see section~\ref{UG:sec:rassi})
\else
(see online manual)
\fi
. We can
start by simply computing the basic interaction properties

%%%To_extract{/doc/samples/problem_based_tutorials/CASSI.excited.acrolein.input}
\begin{inputlisting}
*RASSI calculation on five singlet excited states in acrolein
*File: RASSI.excited.acrolein
*
&GATEWAY
 Title= Acrolein molecule
 coord = acrolein.xyz; basis = STO-3G; group = c1
&SEWARD; &SCF
&RASSCF
 LumOrb
 Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
 CiRoot= 5 5 1
&CASPT2
 Frozen = 4
 MultiState= 5 1 2 3 4 5
                                                                                                                                                                            
>>COPY $Project.JobMix JOB001

&RASSI
 Nr of JobIph
 1 5
 1 2 3 4 5
 EJob
\end{inputlisting}
%%%To_extract

Oscillator strengths for the computed transitions and Einstein coefficients are
compiled at the end of the \program{RASSI} output file. To obtain these values,
however, energy differences have been used which are obtained from the previous
CASSCF calculation. Those energies are not accurate because they do not include
dynamic correlation energy and it is better to substitute them by properly
computed values, such those at the CASPT2 level. This is achieved with the
keyword \keyword{Ejob}. 
\ifmanual
More information is available
in section~\ref{TUT:sec:rassi_thio}.
\fi

Now a more complex case. We want to compute vertical singlet-triplet gaps from
the singlet ground state of acrolein to different, up to five, triplet excited
states. Also, interaction properties are requested. Considering that the spin
multiplicity differs from the ground to the excited states, the spin Hamiltonian
has to be added to our calculations and the \program{RASSI} program takes charge
of that. It is required first, to add in the \program{SEWARD} input the keyword 
\keyword{AMFI}, which introduces the proper integrals required, and to the 
\program{RASSI} input the keyword \keyword{SpinOrbit}. Additionally, as we want
to perform the calculation sequentially and \program{RASSI} will read from
two different wave function calculations, we need to perform specific links
to save the information. The link to the first \program{CASPT2} calculation
will saved in file \file{\$Project.JobMix.S} the data from the \program{CASPT2}
result of the ground state, while the second link before the second \program{CASPT2}
run will do the same for the triplet states. Later, we link these files as
\file{JOB001} and \file{JOB002} to become input files for \program{RASSI}.
In the \program{RASSI} input \keyword{NrofJobIph} will be set to two, meaning
two \file{JobIph} or \file{JobMix} files, the first containing one root (the ground
state) and the second five roots (the triplet states). Finally, we have added 
\keyword{EJob}, which will read the CASPT2 (or MS-CASPT2) energies from the
\file{JobMix} files to be incorporated to the \program{RASSI} results.
The magnitude of properties computed with spin-orbit coupling (SOC) depends
strongly on the energy gap, and this has to be computed at the highest possible
level, such as CASPT2.

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.S-T_gap.acrolein.input}
\begin{inputlisting}
*CASPT2/RASSI calculation on singlet-triplet gaps in acrolein
*File: CASPT2.S-T_gap.acrolein
*
&GATEWAY
 Title= Acrolein molecule
 coord = acrolein.xyz; basis = STO-3G; group= c1
&SEWARD 
 AMFI
&SCF
&RASSCF
 Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
 CiRoot= 1 1 1
&CASPT2
 Frozen= 4
 MultiState= 1 1
>>COPY $Project.JobMix JOB001
&RASSCF
 LumOrb
 Spin= 3; Nactel= 6 0 0; Inactive= 12; Ras2= 5
 CiRoot= 5 5 1
&CASPT2
 Frozen= 4
 MultiState= 5 1 2 3 4 5
>>COPY $Project.JobMix JOB002
&RASSI 
 Nr of JobIph= 2 1 5; 1; 1 2 3 4 5
 Spin
 EJob
\end{inputlisting}
%%%To_extract

As here with keyword \keyword{AMFI}, 
when using command \command{Coord} to build a \program{SEWARD} input
and we want to introduce other keywords, it is enough if we place them
after the line corresponding to \command{Coord}.
Observe that the nature of the triplet states obtained is in sequence one
$n\pi^*$, two $\pi\pi^*$, and two $n\pi^*$. The \program{RASSI} output is 
somewhat complex to analyze, but it makes tables summarizing oscillator 
strengths and Einstein coefficients, if those are the magnitudes of interest. 
Notice that a table is first done with the spin-free states, while the final 
table include the spin-orbit coupled eigenstates (in the CASPT2 energy order 
here), in which each former triplet state has three components.

In many cases working with symmetry will help us to perform calculations
in quantum chemistry. As it is a more complex and delicate problem we direct
the reader to the examples section in this manual. However, we include here
two inputs that can help the beginners. They are based on trans-1,3-butadiene,
a molecule with a C$_{2h}$ ground state. If we run the next input, the
\program{SEWARD} and \program{SCF} outputs will help us to understand how
orbitals are classified by symmetry, whereas reading the \program{RASSCF} output
the structure of the active space and states will be clarified.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.excited.tButadiene.1Ag.input}
\begin{inputlisting}
*CASSCF SA calculation on 1Ag excited states in tButadiene
*File: CASSCF.excited.tButadiene.1Ag
*
&SEWARD 
  Title= t-Butadiene molecule
  Symmetry= Z XYZ
Basis set
C.STO-3G...
C1   -3.2886930 -1.1650250 0.0000000  Bohr
C2   -0.7508076 -1.1650250 0.0000000  Bohr
End of basis
Basis set
H.STO-3G...
H1   -4.3067080  0.6343050 0.0000000  Bohr
H2   -4.3067080 -2.9643550 0.0000000  Bohr
H3    0.2672040 -2.9643550 0.0000000  Bohr
End of basis
                                                                                                                                                                            
&SCF 
                                                                                                                                                                            
&RASSCF 
 LumOrb
 Title= tButadiene molecule (1Ag states). Symmetry order (ag bg bu au)
 Spin= 1; Symmetry= 1; Nactel= 4 0 0; Inactive= 7 0 6 0; Ras2= 0 2 0 2
 CiRoot= 4 4 1

&GRID_IT 
 All
\end{inputlisting}
%%%To_extract

Using the next input will give information about states of a different symmetry.
Just run it as a simple exercise.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.excited.tButadiene.1Bu.input}
\begin{inputlisting}
*CASSCF SA calculation on 1Bu excited states in tButadiene
*File: CASSCF.excited.tButadiene.1Bu
*
&SEWARD 
 Title= t-Butadiene molecule
 Symmetry= Z XYZ
Basis set
C.STO-3G...
C1   -3.2886930 -1.1650250 0.0000000  Bohr
C2   -0.7508076 -1.1650250 0.0000000  Bohr
End of basis
Basis set
H.STO-3G...
H1   -4.3067080  0.6343050 0.0000000  Bohr
H2   -4.3067080 -2.9643550 0.0000000  Bohr
H3    0.2672040 -2.9643550 0.0000000  Bohr
End of basis
                                                                                                                                                                            
&SCF 
                                                                                                                                                                            
&RASSCF 
 FileOrb= $Project.ScfOrb
 Title= tButadiene molecule (1Bu states). Symmetry order (ag bg bu au)
 Spin= 1; Symmetry= 1; Nactel= 4 0 0; Inactive= 7 0 6 0
 Ras2= 0 2 0 2
 CiRoot= 4 4 1
>COPY $Project.RasOrb $Project.1Ag.RasOrb
>COPY $Project.JobIph JOB001

&GRID_IT 
 Name= $Project.1Ag.lus
 All

&RASSCF 
 FileOrb= $Project.ScfOrb
 Title= tButadiene molecule (1Bu states). Symmetry order (ag bg bu au)
 Spin= 1; Symmetry= 3; Nactel= 4 0 0; Inactive= 7 0 6 0; Ras2= 0 2 0 2
 CiRoot= 2 2 1
>COPY $Project.RasOrb $Project.1Bu.RasOrb
>COPY $Project.JobIph JOB002
                                                                                                                                                                            
&GRID_IT 
 Name= $Project.1Bu.lus
 All

&RASSI 
 NrofJobIph= 2 4 2; 1 2 3 4; 1 2
\end{inputlisting}
%%%To_extract

Structure optimizations can be also performed at the CASSCF, RASSCF or CASPT2
levels. Here we shall optimize the second singlet state in the first (here the
only) symmetry for acrolein at the SA-CASSCF level. It is strongly recommended
to use the State-Average option and avoid single state CASSCF calculations for
excited states. Those states are non-orthogonal with the ground state and
are typically heavily contaminated. The usual set of input commands will be
prepared, with few changes. In the \program{RASSCF} input two states will
be simultaneously computed with equal weight (\keyword{CiRoot} 2 2 1), but,
in order to get accurate gradients for a specific root (not an averaged one),
we have to add \keyword{Rlxroot} and set it to two, which is, among the
computed roots, that we want to optimize. The proper density matrix will be
stored. The \program{MCLR} program optimizes, using a perturbative approach,
the orbitals for the specific root (instead of using averaged orbitals), but
the program is called automatically and no input is needed.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.excited_state_optimization.acrolein.input}
\begin{inputlisting}
*CASSCF excited state optimization in acrolein
*File: CASSCF.excited_state_optimization.acrolein
*
 &GATEWAY
Title= acrolein minimum optimization in excited state 2
Basis set
O.STO-3G...2s1p.
O1       1.608542      -0.142162       3.240198 Angstrom
End of basis
Basis set
C.STO-3G...2s1p.
C1      -0.207776       0.181327      -0.039908 Angstrom
C2       0.089162       0.020199       1.386933 Angstrom
C3       1.314188       0.048017       1.889302 Angstrom
End of basis
Basis set
H.STO-3G...1s.
H1       2.208371       0.215888       1.291927 Angstrom
H2      -0.746966      -0.173522       2.046958 Angstrom
H3      -1.234947       0.213968      -0.371097 Angstrom
H4       0.557285       0.525450      -0.720314 Angstrom
End of basis
>>> Do while

 &SEWARD  
                                                                                                                                                                            
>>> If ( Iter = 1 ) <<<
                                                                                                                                                                            
 &SCF 
Title= acrolein minimum optimization
                                                                                                                                                                            
>>> EndIf <<<

 &RASSCF 
LumOrb
Title= acrolein
Spin= 1; nActEl= 4 0 0; Inactive= 13; Ras2= 4
CiRoot= 2 2 1
Rlxroot= 2
                                                                                                                                                                            
 &SLAPAF 
                                                                                                                                                                            
>>> EndDo
\end{inputlisting}
%%%To_extract

In case of performing a \program{CASPT2} optimization for an excited
state, still the SA-CASSCF approach can be used to generate the reference
wave function, but keyword \keyword{Rlxroot} and the use of the \program{MCLR} program
are not necessary, because \program{CASPT2} takes care of selecting
the proper root (the last one).

A very useful tool recently included in \molcas\ is the possibility to
compute minimum energy paths (MEP), representing steepest descendant minimum
energy reaction paths which are built through a series of geometry optimizations, 
each requiring the minimization of the potential energy on a hyperspherical
cross section of the PES centered on a given reference geometry and characterized 
by a predefined radius. One usually starts the calculation from a high energy reference
geometry, which may correspond to the Franck-Condon (FC) structure on an excited-state PES 
or to a transition structure (TS). Once the first lower energy optimized structure is
converged, this is taken as the new hypersphere center, and the procedure is iterated 
until the bottom of the energy surface is reached. Notice that in the TS case a pair of
steepest descent paths, connecting the TS to the reactant and product structures 
(following the forward and reverse orientation of the direction defined by the transition 
vector) provides the minimum energy path (MEP) for the reaction. As mass-weighted 
coordinates are used by default, the MEP coordinate corresponds to the so-called Intrinsic 
Reaction Coordinates (IRC). We shall compute here the MEP from the FC structure of acrolein
along the PES of the second root in energy at the CASSCF level. It is important to remember 
that the CASSCF order may not be accurate and the states may reverse orders at higher
levels such as CASPT2.

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.mep_excited_state.acrolein.input}
\begin{inputlisting}
*CASSCF excited state mep points in acrolein
*File: CASSCF.mep_excited_state.acrolein
*
 &GATEWAY
Title = acrolein mep calculation root 2
Basis set
O.STO-3G...2s1p.
 O1    1.367073     0.000000     3.083333 Angstrom
End of basis
Basis set
C.STO-3G...2s1p.
 C1    0.000000     0.000000     0.000000 Angstrom
 C2    0.000000     0.000000     1.350000 Angstrom
 C3    1.367073     0.000000     1.833333 Angstrom
End of basis
Basis set
H.STO-3G...1s.
 H1    2.051552     0.000000     0.986333 Angstrom
 H2   -0.684479     0.000000     2.197000 Angstrom
 H3   -1.026719     0.000000    -0.363000 Angstrom
 H4    0.513360     0.889165    -0.363000 Angstrom
End of basis

>>> EXPORT MOLCAS_MAXITER=300
>>> Do while

 &SEWARD
>>> If ( Iter = 1 ) <<<
 &SCF 
>>> EndIf <<<

 &RASSCF 
   Title="acrolein mep calculation root 2"; Spin=1
   nActEl=4 0 0; Inactive=13; Ras2=4; CiRoot=2 2 1; Rlxroot=2
 &SLAPAF 
   MEP-search
   MEPStep=0.1

>>> EndDo
\end{inputlisting}
%%%To_extract

As observed, to prepare the input for the MEP is simple, just add the keyword \keyword{MEP-search}
and specify a step size with \keyword{MEPStep}, and the remaining structure equals that of a geometry optimization.
The calculations are time consuming, because each point of the
MEP (four plus the initial one obtained here) is computed through a specific optimization.
A file named \file{\$Project.mep.molden} (read by \program{MOLDEN} )
will be generated in \$WorkDir containing only those points belonging to the MEP.

We shall now show how to perform geometry optimizations under nongeometrical
restrictions, in particular, how to compute hypersurface crossings, which are key structures
in the photophysics of molecules. We shall get those points as minimum energy crossing points in
which the energy of the highest of the two states considered is minimized under the restriction
that the energy difference with the lowest state should equal certain value (typically zero).
Such point can be named a minimum energy crossing point (MECP). If a further restriction is
imposed, like the distance to a specific geometry, and several MECP as computed at varying distances,
it is possible to obtain a crossing seam of points where the energy between the two states is
degenerated. Those degeneracy points are funnels with the highest probability for the energy
to hop between the surfaces in internal conversion or intersystem crossing photophysical processes.
There are different possibilities. A crossing between states of the same spin
multiplicity and spatial symmetry is named a conical intersection. Elements like the nonadiabatic
coupling terms are required to obtain them strictly, and they are not computed presently
by \molcas. If the crossing occurs between states of the same
spin multiplicity and different spatial symmetry or between states of different spin multiplicity,
the crossing is an hyperplane and its only requirement is the energetic degeneracy and the
proper energy minimization.

Here we include an example with the crossing between the lowest singlet (ground) and triplet 
states of acrolein. Notice that two different states are computed, first by using 
\program{RASSCF} to get the wave function and then \program{ALASKA} to get the gradients
of the energy. Nothing new on that, just the information needed in any geometry optimizations.
The \program{GATEWAY} input requires to add as constraint an energy
difference between both states equal to zero. A specific instruction is required after 
calculating the first state. We have to copy the communication file \file{RUNFILE}
(at that point contains the information about the first state) to \file{RUNFILE2}
to provide later \program{SLAPAF} with proper information about both states:

%%%To_extract{/doc/samples/problem_based_tutorials/CASSCF.S-T_crossing.acrolein.input}
\begin{inputlisting}
*CASSCF singlet-triplet crossing in acrolein
*File: CASSCF.S-T_crossing.acrolein
*
 &GATEWAY
Title= Acrolein molecule
Basis set
O.sto-3g....
 O1             1.5686705444       -0.1354553340        3.1977912036  Angstrom
End of basis
Basis set
C.sto-3g....
 C1            -0.1641585340        0.2420235062       -0.0459895824  Angstrom
 C2             0.1137722023       -0.1389623714        1.3481527296  Angstrom
 C3             1.3218729238        0.1965728073        1.9959513294  Angstrom
End of basis
Basis set
H.sto-3g....
 H1             2.0526602523        0.7568282320        1.4351034056  Angstrom
 H2            -0.6138178851       -0.6941171027        1.9113821810  Angstrom
 H3            -0.8171509745        1.0643342316       -0.2648232855  Angstrom
 H4             0.1260134708       -0.4020589690       -0.8535699812  Angstrom
End of basis
Constraints
   a = Ediff
  Value
   a = 0.000
End of Constraints

>>> Do while
                                                                                                                                                                            
 &SEWARD  
                                                                                                                                                                            
>>> IF ( ITER = 1 ) <<<
 &SCF 
>>> ENDIF <<<
                                                                                                                                                                            
 &RASSCF 
   LumOrb
   Spin= 1; Nactel= 4 0 0; Inactive= 13; Ras2= 4
   CiRoot= 1 1; 1
 &ALASKA
>>COPY $WorkDir/$Project.RunFile $WorkDir/RUNFILE2
                                                                                                                                                                            
 &RASSCF 
   LumOrb
   Spin= 3; Nactel= 4 0 0; Inactive= 13; Ras2= 4
   CiRoot= 1 1; 1
 &ALASKA 
 &SLAPAF 
                                                                                                                                                                            
>>> EndDo
\end{inputlisting}
%%%To_extract

%$
Solvent effects can be also applied to excited states, but first the reaction
field in the ground (initial) state has to be computed. This is because solvation in
electronic excited states is a non equilibrium situation in with the electronic
polarization effects (fast part of the reaction field) have to treated apart
(they supposedly change during the excitation process) from the orientational
(slow part) effects. The slow fraction of the reaction field is maintained from
the initial state and therefore a previous calculation is required. 
From the practical point of view the input is simple as illustrated in the next 
example. First, the proper reaction-field
input is included in \program{SEWARD}, then a \program{RASSCF} and \program{CASPT2}
run of the ground state, with keyword \keyword{RFPErt} in \program{CASPT2},
and after that another SA-CASSCF calculation of five roots to get the wave function
of the excited states. Keyword \keyword{NONEequilibrium} tells the program to extract
the slow part of the reaction field from the previous calculation of the ground
state (specifically from the \file{JOBOLD} file, which may be stored for other
calculations) while the fast part is freshly computed. Also, as it is a SA-CASSCF
calculation (if not, this is not required) keyword \keyword{RFRoot} is introduced
to specify for which of the computed roots the reaction field is generated. We have
selected here the fifth root because it has a very large dipole moment, which is
also very different from the ground state dipole moment. If you compare the excitation
energy obtained for the isolated and the solvated system, a the large red shift is 
obtained in the later.

%%%To_extract{/doc/samples/problem_based_tutorials/CASPT2.excited_solvent.acrolein.input}
\begin{inputlisting}
*CASPT2 excited state in water for acrolein
*File: CASPT2.excited_solvent.acrolein
*
&GATEWAY  
  Title= Acrolein molecule
  coord = acrolein.xyz; basis = STO-3G; group= c1
  RF-input
   PCM-model; solvent= water
  End of RF-input
&SEWARD  
&RASSCF 
  Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
  CiRoot= 1 1 1
&CASPT2 
  Multistate= 1 1
  RFPert
&RASSCF 
  Spin= 1; Nactel= 6 0 0; Inactive= 12; Ras2= 5
  CiRoot= 5 5 1
  RFRoot= 5
  NONEquilibrium
&CASPT2 
  Multistate= 1 5
  RFPert
\end{inputlisting}
%%%To_extract

A number of simple examples as how to proceed with the most frequent 
quantum chemical problems computed with \molcas\ have been given above. Certainly there are many more
possibilities in \molcas\ \molcasversion\, such as calculation of 3D band
systems in solids at a semiempirical level, obtaining valence-bond structures,
the use of QM/MM methods in combination with a external MM code, the introduction
of external homogeneous or non homogeneous perturbations, generation of atomic
basis sets, application of different localization schemes, analysis of first
order polarizabilities, calculation of vibrational intensities, analysis, generation,
and fitting of potentials, computation of vibro-rotational spectra for diatomic
molecules, introduction of relativistic effects, etc. All those aspects are
explained in the manual and are much more specific. Next section~\ref{TUT:sec:pg-based-tut}
details the basic structure of the inputs, program by program, while easy examples
can also be found. Later, another chapter includes a number of extremely detailed
examples with more elaborated quantum chemical examples, in which also scientific
comments are included. Examples include calculations on high symmetry molecules,
geometry optimizations and Hessians, computing reaction paths, high quality wave
functions, excited states, solvent models, and computation of relativistic effects.

