|
MOLCAS manual: Next: 5.3 Computing excited states Up: 5. Problem Based Tutorials Previous: 5.1 Electronic Energy at Fixed
5.2 Optimizing geometriesIt 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 8.1 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
The general input commands Do while and 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.
(see section All examples previously discussed, use COORD keyword, but it also possible to use native format, where symmetry unique atoms are specified (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:
The sequence of programs employed includes GATEWAY which is external to the loop, followed by SEWARD, SCF, and SLAPAF. SEWARD computes the integrals, SCF program computes the RHF energy and wave function, and 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 $Project.geo.molden will be generated in $WorkDir containing all the geometric steps contained in the optimization process. MOLDEN or 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 H2O+ cation:
The same procedure can be followed if we pretend to perform a DFT geometry optimization:
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 MCKINLEY computes second derivatives of a predefined (SCF or CASSCF) wave function, while MCLR performs the vibrational and statistical analyses. MOLCAS simply requires input for the MCKINLEY program to perform the entire calculation by using keywords Perturbation and Hessian, while 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.
Note that MCKINLEY input above is placed after 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 SEWARD input is present. To avoid that, the only solution is to remove the communication 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 $Project.freq.molden is generated containing the vibrational frequencies and modes, which can be visualized by 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 120o
during the optimization.
This means that only the O-H bond distances be optimized in this partial minimization.
The restriction is indicated
in in GATEWAY
by invoking the keyword Constraints and ending with the 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.
(see nomenclature in
section Note that the RASSCF program requires initial trial orbitals, and those which are automatically generated by SEWARD are used. The resulting CASSCF wave function includes all valence orbitals and electrons.
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 Å), while the
bond angle, having an initial of 81
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 112o. 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 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. An illustration of the input required for transition state optimization for water at the DFT level is given below:
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 C2v symmetry. A transition state has been found with a linear HOH angle of 180o. 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 FindTS in SLAPAF. FindTS must 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, 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 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 NUMERICAL_GRADIENT, which will take care of performing numerical gradients. From the user point of view the only requirement 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-12+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 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, 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 NUMERICAL_GRADIENT
has been reduced to five instead of seven, because many of the deformations
are redundant within the C2v symmetry. Also, symmetry considerations are
important when defining geometrical restrictions
(see sections
    
Next: 5.3 Computing excited states Up: 5. Problem Based Tutorials Previous: 5.1 Electronic Energy at Fixed |
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.
![[*]](crossref.png)
, is optimized.
