|
MOLCAS manual: Next: 4.2 Program-based tutorials Up: 4. Tutorials Previous: 4. Tutorials
4.1 Problem-based tutorials
4.1.1 Computing the electronic energy at fixed nuclear geometryThe MOLCAS 7.4 suite of quantum chemical programs is modular in design. The 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 thefunctions 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. Our initial problem will be to compute an electronic energy at a fixed geometry of the nuclei, and this will be done using different methods and thus different program modules. First, the proper MOLCAS environment has to be set up. The following variables must be properly defined, for instance:
If not defined, MOLCAS gives them the default values. The MOLCAS variable will be set up as the latest implemented version of the code. Once in the MOLCAS home directory the variable is set up directly. Project and WorkDir have the default values None and $PWD, respectively. It is very important that the molcas driver, called by command molcas and built during the installation of the code, is included in the $PATH. We prepare our first run: the calculation of the SCF energy of the methane (CH4) molecule. Three programs should be used: GATEWAY to specify information about the system, SEWARD to compute and store the one- and two-electron integrals, and SCF to obtain the Hartree-Fock SCF wave function and energy. Three inputs should be prepared. The GATEWAY input will contain the nuclear geometry in cartesian coordinates and the label for the one-electron basis set, here, for instance, an STO-3G minimal basis set. We shall use here the keyword coord, which allows to automatically create a GATEWAY input from a standard file containing the cartesian coordinates in Angstrom (as generated by many programs like MOLCAS GV or MOLDEN). No symmetry is being considered so far, therefore we introduce the keyword nosym to force the program not to look for any symmetry in the molecule. In this case the input for SEWARD is void. In a closed-shell case like this one the SCF input can be left empty. All the input files used here can be found at $MOLCAS/doc/examples/tutorials, like the file SCF.energy.CH4 described below.
where the content of the CH4.xyz file is:
Now, to run MOLCAS we simply type:
and the main output will be stored in file SCF.energy.CH4.out (or in the second case SCF.energy.CH4.log: and the default error file in SCF.energy.CH4.err. The most relevant information will be contained in the output file, where we shall analyze the SEWARD 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 4.2.2. Next, comes the output of program SCF with information of the electronic energy, wave function, and the Hartree-Fock (HF) molecular orbitals (see section 4.2.3). Files containing intermediate information, integrals, orbitals, etc, will be kept in the $WorkDir directory for further use. For instance files $Project.OneInt and $Project.OrdInt contains the one- and two-electron integrals stored in binary format. File $Project.ScfOrb keeps the HF molecular orbitals in ASCII format. Finally, $Project.RunFile is a communication file between programs. All those files can be later used for more advanced calculations in order to avoid repeating the already run calculations.
There are some graphical utilities that can be used for the analysis of the
results. By default, MOLCAS generates files which can be read with the
MOLDEN program. In particular we can find in $WorkDir the file
CH4.scf.molden, where information about molecular geometry, molecular
orbitals (use Density Mode in MOLDEN), etc, can be found (See also
section
while option gv -c should be used if reading from a SEWARD input. Once in the screen the molecule can be manipulated in many ways, fo example, changing its geometry, adding or deleting atoms, etc. You can also use it to build new molecules. In order to obtain the information to display molecular orbitals and densities it is necessary to run a MOLCAS program named GRID_IT
and run it:
GRID_IT can also be run separately if an orbital file is specified in the input and the $WorkDir directory is available. In the $WorkDir and $PWD directories a new file is generated, CH4.grid which contains the information required by the GRID_IT input. By typing:
a window will be opened displaying the molecule and its charge density. By proper
selection with the mouse buttons, shape and size of several molecular orbitals
can be selected. Typing molcas gv -? one can get a short help for the
use of the program. More information can be found in sections As an alternative to run the projects we include here a short script to be placed in the directory $MOLCAS/doc/samples/problem_based_tutorials as the file project.sh. Just run project.sh $Project, where $Project is the MOLCAS input and you will get output and error files and a $WorkDir directory named $Project.work.
In order to run a Kohn-Sham density functional calculation, MOLCAS uses the same SCF module, therefore the only change needed is to specify in the SCF input the DFT option and the required functional, for instance here, B3LYP:
and the same graphical files can be found in $WorkDir and $PWD. Our next step is to get the second-order Møller-Plesset perturbation (MP2) energy for methane at the same molecular geometry and using the same one-electron basis set. Program MBPT2 has to be invoked this time. We can take advantage of having previously computed the proper integrals with SEWARD and the reference closed-shell HF wave function with the SCF program. In such case we may keep the same definitions as before and simply prepare a file containing the MBPT2 input and run it using the molcas command. The proper intermediate file will be already in $WorkDir. If, on the contrary, one has to start from scratch, all required inputs should be placed sequentially, for instance, in the MP2.energy.CH4 file. If we decide to start the project from the beginning is convenient to remove the $WorkDir directory, unless we are certain about the files we are going to find there.
On top of the HF wave function, an MP2 calculation has been performed with the deepest orbital, carbon 1s, of the molecule frozen. Information about the output of the MBPT2 program can be found on section 4.2.13. The SCF program works by default with closed-shell systems with an even number of electrons at the Restricted Hartee-Fock (RHF) level. We may instead want to use the Unrestricted Hartree Fock (UHF) method, by invoking the keyword UHF. This is possible for both even and odd electron systems. For instance in a system with odd number of electrons like the CH3 radical
the input to run an open-shell UHF calculation will be simply
If we want to force a charged system, this must be indicated in the SCF input, for instance to compute the cation of the CH4 molecule at the UHF level:
The Kohn-Sham DFT calculation can be also run on the basis of the UHF algorithm
In the UHF (and UHF/DFT) method it is also possible to force particular
In the output the final occupations are six Different sets of methods use other MOLCAS modules. For instance, to perform a Complete Active Space (CAS) SCF calculation, the RASSCF program has to be used. This module requires starting trial orbitals, which can be obtained from a previous SCF calculation or automatically by SEWARD, which provides trial orbitals by using a model Fock operator. Recommended keywords are Nactel, defining the total number of active electrons, holes in Ras1, and particles in Ras3, respectively (the two last ones only for RASSCF-type calculations), Inactive, including the number of orbitals inactive (occupation always two) in the CASSCF reference, and Ras2, defining the number of active orbitals. By default the obtained wave function will correspond to the lowest state of the symmetry with spin multiplicity one. Most of the input can be avoided if one has prepared and linked an INPORB file with the different orbital types defined (for example with the programs GV).
In this case we are computing the lowest singlet state, the ground state here, considering that this is a closed-shell situation, with an active space of eight electrons in eight orbitals. The lowest orbital of the molecule (carbon 1s) has been left inactive. This is a CASSCF case in which all the valence orbitals and electrons (carbon 2s2p plus four hydrogen 1s) have been included in the active space. Using the CASSCF wave function as a reference, we can perform a second-order perturbative, CASPT2, correction to the electronic energy by employing the CASPT2 program. If we maintain in the $WorkDir directory the files generated by the previous run (in particular the following files: integral files CH4.OneInt,CH4.OrdInt; CASSCF wave function information file CH4.JobIph, and communication file CH4.RunFile), it will not be necessary to run again programs SEWARD, and RASSCF. In that case it will be enough to prepare an input for the CASPT2 program and run it. Here, however, we include the full input file:
As an alternative to generate trial orbitals for the RASSCF program, we have added here an input for the SCF program. In most of cases the Hartree-Fock orbitals will be a better choice as starting orbitals. In that case, the RASSCF input has to include keyword LumOrb to read from any external source of orbitals other than those generated by the SEWARD program. We have also changed the input of the RASSCF program. Now the number of active orbitals is six, as well as the active electrons. We have excluded from the space the two lowest orbitals (Inactive 2) and also another orbital goes to the secondary space. If we had stayed with the previous (8,8) full valence space, considering that we are using a minimal basis set, the CASPT2 program would not have more electronic correlation energy to include. Regarding the input for the CASPT2 program we have frozen the lowest orbital (carbon 1s). We may change the charge and multiplicity of our wave function by computing the CH4+ cation with the same methods. It is the RASSCF program which defines the character of the problem by specifying number of electrons, and spin and spatial symmetry. We have here one electron less and a doublet multiplicity.
No further modification is needed in the CASPT2 input:
A somewhat more sophisticated calculation can be performed at the Restricted Active Space (RAS) SCF level. In such case the level of excitation in the CI expansion can be controlled by restricting the number of holes and particles present in certain orbitals.
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, either coming from Ras1 or Ras2. Within Ras2 all type of orbital occupations are allowed. The RASSCF wave functions can, in principle, be used as reference for multiconfigurational perturbation theory (RASPT2) but this approach is still in the experimental stage. MOLCAS also has the possibility of computing electronic energies at different CI levels by using the MRCI program. We have prepared here a Singles and Doubles CI (SDCI) calculation on the CH4 molecule. To set up the calculations, programs MOTRA, which transforms the integrals to molecular basis, and GUGA, which computes the coupling coefficients, have to be run before the MRCI program. In MOTRA we have specified the reference orbitals employed (those from an HF SCF calculation) and the frozen orbitals. In GUGA we describe the reference for the CI, with the number of correlated electrons, the spatial and spin symmetry, the inactive orbitals (always occupation two in the reference space), and the type of CI expansion.
If we want to use reference orbitals from a previous CASSCF calculation, the RASSCF program will have to be run before the MOTRA module. Also, if the spatial or spin symmetry are changed for the CI calculation, the modifications will be introduced in the GUGA input. Many alternatives are possible in order to perform an MRCI calculation like next one, in which the reference space to perform the CI is multiconfigurational:
The MRCI program allows also to get electronic energies using the ACPF method. Another MOLCAS program, CPF, offers the possibility to use the CPF, MCPF, and ACPF methods with a single reference function. The required input is quite similar than that for the MRCI program:
Finally, MOLCAS can also perform closed- and open-shell coupled cluster calculations at the CCSD and CCSD(T) levels. They are controlled by the CCSDT program, whose main requirement is that the reference function has to be generated with RASSCF. The following input is required to obtain the CCSD(T) energy for the CH4 molecule:
As it is a closed-shell calculation, the RASSCF input is prepared to compute a simple RHF wave function (zero active electrons and orbitals), with keywords OutOrbitals Canonical. The MOTRA has to include the keyword JobIph to extract the wave function information from file JOBIPH (automatically generated by RASSCF). Finally, CCT in program CCSDT leads to the calculation of the CCSD(T) energy using the default algorithms. The CCSDT program in MOLCAS is specially suited to compute open-shell cases. The input required to obtain the electronic energy of the CH4+ cation with the CCSD(T) method is:
where it is the RASSCF program which generated the proper Restricted Open-Shell (RO) HF reference. Different levels of spin adaptation are available. We may want to introduce the effect of solvents into our calculation. MOLCAS includes two models: Kirkwood and PCM. To add solvent effect to a ground state at the HF, DFT, or CASSCF levels simply requires to include the RF-input within SEWARD, while the reaction field is computed self-consistently.
Other programs such as CASPT2, RASSI, and MOTRA require that the reaction field is included as a perturbation with keyword RFPErturbation. In the next example the correction is added at both the CASSCF and CASPT2 levels.
Notice that the tesserae of the average area in the PCM model (keyword
Aare has been changed to the value required for acetone, while the
default is 0.4 Å2 for water (see section
4.1.2 Optimizing geometries: minima, transition states, crossings, and minimum energy pathsWe may want now to explore potential energy surfaces (PES) and optimize the molecular geometry for specific points in the PES. Different cases can be accomplished, like getting 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 in which the energy becomes degenerate, or mapping the minimum (steepest descendent) energy path (MEP) from a specific point if the PES downward to a minimum energy situation. All these searches can be performed fully optimizing all the degrees of freedom of the system or introducing certain restrictions. MOLCAS 7.4 can perform geometry optimizations at the SCF (RHF and UHF), DFT (RHF and UHF based), CASSCF (CASSCF and RASSCF), and CASPT2 and other correlated levels, although in the last case only numerical gradients are available, where in the other cases more efficient analytical gradients are used.
Geometry optimizations require many loops, in which the electronic energy is estimated at a specific
level of calculation and the gradients of the energy (first derivatives) with respect to the geometric
degrees of freedom (DOF) are computed. With this information at hand the program has to decide if we are
already at the final required geometry (for instance if gradient is zero for all
DOF we may have reached a minimum in the PES) or if we must change the geometry
to try to get a better solution in a new iteration. The input file should
therefore be built in a way that allows to loop over the different programs. The
general input commands Do while and Enddo control the loop
and the input is be inserted within them. We may want to add also instructions
about the number of maximum iterations allowed and the type of output required
(see section The selected example is the geometry optimization of the water molecule at the SCF RHF level of calculation:
The sequence of programs employed is GATEWAY (out of the loop), and then SEWARD, SCF, and SLAPAF. SEWARD is computes the integrals SCF program computes the RHF energy, and wave function. SLAPAF will control the calculation of gradientns and estimate if the calculation has already finished or need to proceed to a new nuclear geometry for the next iteration. Automatically a file named $Project.geo.molden will be generated in $WorkDir containing all the geometric steps contained in the optimization process. MOLDEN can read such file. To move to another reference wave function means simply to change the initial programs. For instance, we can perform an UHF calculation of the H2O+ cation:
The same procedure can be followed if we pretend to perform a DFT geometry optimization:
Once we have found an energy minimum based on the calculation of gradients, we may ensure that this is really a minimum energy point, what can be only accomplished by computing second derivatives of the energy, that is, the Hessian. MOLCAS can compute analytical Hessians for SCF and single state CASSCF wave functions. For some other methods one can use numerical procedures to compute the Hessian. Once the Hessian is computed we have the vibrational frequencies and we can apply statistical mechanics to obtain thermodynamic properties. If we are in a true energy minimum, all frequencies will be real (instead of imaginary). Program MCKINLEY computes second derivatives of a predefined (SCF or CASSCF) wave function, while MCLR performs the vibrational and statistical analyses. MOLCAS simply requires an input for the MCKINLEY program to do all the job by using keywords Perturbation Hessian, while program MCLR will be called automatically and no input is required. We include here the full set of calculations, first, a geometry optimization, followed by the calculation of a Hessian.
Notice that MCKINLEY input is placed after EndDo, therefore, out of the looping scheme. Once the geometry optimization has finished, the Hessian will be computed at the final obtained geometry and computation level. In general, any calculation performed using a $WorkDir directory where a previous geometry optimization of any type has taken place will use as geometry the last one obtained in the optimization, even if a SEWARD input is present. To avoid that, the only solution is to remove the communication file RUNFILE where the geometry is stored. Notice also that the frequencies are computed in a cartesian basis, and that three translational and three rotational frequencies (they 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 in its minimum energy structure we expect (3N-6) three real vibrational frequencies. By default in $WorkDir a file $Project.freq.molden is generated containing information about the vibrational frequencies, and modes, which can be visualized by MOLDEN.
We can introduce here a new level of theory, CASSCF, although this is especially suited for
geometry optimizations of excited states (we shall consider excited states in the next chapter).
A geometry optimization is performed, and to show more possibilities this time
we shall impose the restriction that the HOH angle in water should be 120o
at the final structure. That means that only the bond distances will be really
optimized in this partial minimization. The restriction is indicated
in SLAPAF opening the keyword Constraints, which ends with another keyword,
End of Constraints. Inside we assign the name of some variables (one here)
to the geometrical parameters we want to restrict (see nomenclature in
section Observe also that as the RASSCF program requires initial trial orbitals, we use those which were automatically generated from SEWARD. The defined CASSCF wave function includes here all valence orbitals and electrons.
Other more flexible ways to impose geometric restrictions is to decide which internal
coordinates should remain fixed and which should change. In the next example we force
the bond lengths to remain fixed at their initial distance (here 0.91 Å), while the
bond angle (initially 81
It the final output the bond lengths remain as initially while the angle is optimized to reach 112o. Our next step is to compute a transition state, a structure connecting different regions of the potential energy hypersurface which is a maximum just for 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 keyword TS to the SLAPAF input. Keyword PRFC is suggested in order to verify the nature of the transition structure. Searching for transition states is, however, not an easy task. Here we illustrate it for water at the DFT level:
Notice that we have used the initial coordinates in Bohr, just to show that this is the default case. The optimal geometry for ground state water has C2v symmetry. A transition state has been found with a linear HOH angle of 180o. In many cases, as for example along the energy path for a chemical reaction, we may have a clue about the structure of the transition state, which typically represents an intermediate conformation between reactives and products. In that case we may try to orient the program leading it to the guessed solution. Keyword FindTS in SLAPAF has this goal. It has to be accompanied with a definition of constrained geometric definitions. 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 we show an example in which the initial geometry of water is clearly bent, and we impose the trial restriction that the angle for the transition state should be near 180o. The final transition state will, however, be obtained without any type of geometrical restriction.
The CASPT2 geometry optimizations are somewhat different because ALASKA is not suited to compute CASPT2 analytical gradients. Therefore the ALASKA program is automatically substituted by program CASPT2_GRADIENT, which will take care of performing numerical gradients. From the user pointview the only requirent is to place the CASPT2 input after the 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+1 CASPT2 calculations are performed each pass of the optimization, where N is the number of atoms.
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 CASPT2 optimization by restricting the molecule to work in the C2v point group, which, by the way, is the proper symmetry for water in the ground state. The GATEWAY program (as no symmetry has been specified) will identify and work with the highest available point group, C2v. Here the molecule is placed with YZ as the molecular plane. By adding keyword Symmetry containing as elements of symmetry the YZ (symbol X) and YX (symbol Z), the poing 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, a1, b1, b2, and a2, 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 SEWARD and must be checked before proceeding, because from that point the elements of symmetry will be known by their order in SEWARD: a1, b1, b2, and a2, for instance, will be symmetries 1, 2, 3, and 4, respectively. SCF does not require to specify the class of orbitals and it can be used as a learning tool.
Thanks to symmetry restrictions the number of iterations within CASPT2_GRADIENT
has been reduced to nine instead of nineteen, because many of the deformations
are redundant within the C2v symmetry. Also, symmetry considerations are
important when defining geometrical restrictions
(see sections
4.1.3 Computing excited statesThe 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 (see figure above) at the CASSCF level and using a minimal basis set. The standard file with cartesian coordinates is:
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.
We have used as active all the
In 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
GRID_IT. If we use now the GV program for instance with the
command molcas gv $Project.grid -a -1.0 $Project.grid2 the graph
will display the differential charge transfer from the ground to the
first excited state, in which the charge will leave the oxygen lone pair
toward the 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 SEWARD and RASSCF (unnecessary if they were computed previously), adding a CASPT2 input with the keyword MultiState set to 5 1 2 3 4 5. The CASPT2 will perform four consecutive single-state (SS) CASPT2 calculations using the SA-CASSCF roots computed by the 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 FROZen keyword is put here as a reminder. By default the program leaves the core orbitals frozen.
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 RASSI module, a powerful
program which can be used for many purposes (see section
Oscillator strengths for the computed transitions and Einstein coefficients are compiled at the end of the 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 Ejob. More information is available in section 5.5.1. 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 RASSI program takes charge of that. It is required first, to add in the SEWARD input the keyword AMFI, which introduces the proper integrals required, and to the RASSI input the keyword SpinOrbit. Additionally, as we want to perform the calculation sequentially and RASSI will read from two different wave function calculations, we need to perform specific links to save the information. The link to the first CASPT2 calculation will saved in file $Project.JobMix.S the data from the CASPT2 result of the ground state, while the second link before the second CASPT2 run will do the same for the triplet states. Later, we link these files as JOB001 and JOB002 to become input files for RASSI. In the RASSI input NrofJobIph will be set to two, meaning two JobIph or JobMix files, the first containing one root (the ground state) and the second five roots (the triplet states). Finally, we have added EJob, which will read the CASPT2 (or MS-CASPT2) energies from the JobMix files to be incorporated to the 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.
As here with keyword AMFI,
when using command Coord to build a SEWARD input
and we want to introduce other keywords, it is enough if we place them
after the line corresponding to Coord.
Observe that the nature of the triplet states obtained is in sequence one
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 C2h ground state. If we run the next input, the SEWARD and SCF outputs will help us to understand how orbitals are classified by symmetry, whereas reading the RASSCF output the structure of the active space and states will be clarified.
Using the next input will give information about states of a different symmetry. Just run it as a simple exercise.
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 RASSCF input two states will be simultaneously computed with equal weight (CiRoot 2 2 1), but, in order to get accurate gradients for a specific root (not an averaged one), we have to add 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 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.
In case of performing a CASPT2 optimization for an excited state, still the SA-CASSCF approach can be used to generate the reference wave function, but keyword Rlxroot and the use of the MCLR program are not necessary, because CASPT2 takes care of selecting the proper root. A very useful tool recently included in MOLCAS is the possibility to compute minimum energy paths (MEP), representing steepest descendent 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.
As observed, to prepare the input for the MEP is simple, just to add the keyword MEP-search and specify a hypersphere constraint to the SLAPAF input, 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 $Project.mep.molden (read by 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 RASSCF to get the wave function and then ALASKA to get the gradients of the energy. Nothing new on that, just the information needed in any geometry optimizations. The SLAPAF 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 RUNFILE (at that point contains the information about the first state) to RUNFILE2 to provide later SLAPAF with proper information about both states:
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 SEWARD, then a RASSCF and CASPT2 run of the ground state, with keyword RFPErt in CASPT2, and after that another SA-CASSCF calculation of five roots to get the wave function of the excited states. 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 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 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.
A number of simple examples as how to proceed with the most frequent quantum chemical problems computed with MOLCAShave been given above. Certainly there are many more possibilities in MOLCAS 7.4 such as calculation of 3D band systems in solids at a semiempirical level, obaining 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 4.2 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.
    Next: 4.2 Program-based tutorials Up: 4. Tutorials Previous: 4. Tutorials |
MOLCAS=/home/molcas/molcas.7.4\
![[*]](crossref.png)
and 
orbital occupations. There are two ways to do it. First,
with keyword ZSPIn in the SCF program, which represents the
difference between the number of 
) is optimized.
and 
orbitals, two bonding and
two antibonding 
states, and two 
states.
It is convenient to add the GRID_IT input in order to be able to use
the GV 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:
, two 
, and two 
