nextuppreviouscontents
MOLCAS manual:

Next: 10.3 Computing a reaction Up: 10. Examples Previous: 10.1 Computing high symmetry

Subsections



10.2 Geometry optimizations and Hessians.

To optimize a molecular geometry is probably one of the most frequent interests of a quantum chemist [208]. In the present section we examine some examples of obtaining stationary points on the energy surfaces. We will focus in this section in searching of minimal energy points, postponing the discussion on transition states to section [*]. This type of calculations require the computation of molecular gradients, whether using analytical or numerical derivatives. We will also examine how to obtain the full geometrical Hessian for a molecular state, what will provide us with vibrational frequencies within the harmonic approximation and thermodynamic properties by the use of the proper partition functions.

The program ALASKA computes analytical gradients for optimized wave functions. In 7.8 the SCF, DFT, and CASSCF/RASSCF levels of calculation are available. The program ALASKA also computes numerical gradients from CASPT2 and MS-CASPT2 energies. Provided with the first order derivative matrix with respect to the nuclei and an approximate guess of the Hessian matrix, the program SLAPAF is then used to optimize molecular structures. From MOLCAS-5 it is not necessary to explicitly define the set of internal coordinates of the molecule in the SLAPAF input. Instead a redundant coordinates approach is used. If the definition is absent the program builds its own set of parameters based on curvature-weighted non-redundant internal coordinates and displacements [209]. As they depend on the symmetry of the system it might be somewhat difficult in some systems to define them. It is, therefore, strongly recommended to let the program define its own set of non-redundant internal coordinates. In certain situations such as bond dissociations the previous coordinates may not be appropriate and the code directs the user to use instead Cartesian coordinates, for instance.

10.2.1 Ground state optimizations and vibrational analysis

As an example we are going to work with the 1,3-cyclopentadiene molecule. This is a five-carbon system forming a ring which has two conjugated double bonds. Each carbon has one attached hydrogen atom except one which has two. We will use the CASSCF method and take advantage of the symmetry properties of the molecule to compute ground and excited states. To ensure the convergence of the results we will also perform Hessian calculations to compute the force fields at the optimized geometries.

In this section we will combine two types of procedures to perform calculations in MOLCAS. The user may then choose the most convenient for her/his taste. We can use an general script and perform an input-oriented calculation, when all the information relative to the calculation, including links for the files and control of iterations, are inserted in the input file. The other procedure is the classical script-oriented system used in previous examples and typically previous versions of MOLCAS. Let's start by making an input-oriented optimization. A script is still needed to perform the basic definitions, although they can be mostly done within the input file. A suggested form for this general script could be:



#!/bin/sh
export MOLCAS=/home/molcas/molcashome
export MOLCASMEM=64
export Project=Cyclopentadiene1
export HomeDir=/home/somebody/somewhere
export WorkDir=$HomeDir/$Project
[ ! -d $WorkDir ] && mkdir $WorkDir
molcas $HomeDir/$Project.input >$HomeDir/$Project.out 2>$HomeDir/$Project.err
exit

We begin by defining the input for the initial calculation. In simple cases the optimization procedure is very efficient. We are going, however, to design a more complete procedure that may help in more complex situations. It is sometimes useful to start the optimization in a small size basis set and use the obtained approximate Hessian to continue the calculation with larger basis sets. Therefore, we will begin by using the minimal STO-3G basis set to optimize the ground state of 1,3-cyclopentadiene within C2v symmetry.

Figure 10.2: 1,3-cyclopentadiene
\begin{figure}{---------------------------------------------------}
\myincludegraphics{advanced.examples/cyclope}\end{figure}

We will use the following input in an input-oriented calculation. Notice that we have directed the output files sequentially (one per iteration) to the $WorkDir directory by using the Set Output File command, the maximum number of iterations of the subsequent loops, and the starting and end of the loops on each step of the optimization procedure by using the commands Do while and EndDo. It is important than the parameter MaxIter never goes beyond the number of iterations in the SLAPAF input.



>>>  Set  Output  File  <<<
>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

>>LINK  $CurrDir/$Project.ForceConstant.STO-3G  RUNFILE

  &SEWARD  &END
Title
1,3,-cyclopentadiene.  STO-3G  basis  set.
Symmetry
  X  XY
Basis  set
C.STO-3G....
C1  0.000000  0.000000  0.000000  Bohr
C2  0.000000  2.222644  1.774314  Bohr
C3  0.000000  1.384460  4.167793  Bohr
End  of  basis
Basis  set
H.STO-3G....
H1  1.662033  0.000000  -1.245623  Bohr
H2  0.000000  4.167844  1.149778  Bohr
H3  0.000000  2.548637  5.849078  Bohr
End  of  basis
End  of  Input

  &SCF  &END
TITLE
  cyclopentadiene  molecule
OCCUPIED
9  1  6  2
ITERATIONS
40
END  OF  INPUT

  &RASSCF  &END
TITLE
  cyclopentadiene  molecule  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  9  0  6  0
RAS2
  0  2  0  3  <---  All  pi  valence  orbitals  active
ITER
50,25
CIMX
25
LUMORB
END  OF  INPUT

  &ALASKA  &END
End  of  Input

  &SLAPAF  &END
Iterations
80
Thrs
0.5D-06  1.0D-03
End  of  Input

>>>  EndDo  <<<

A link to the RUNFILE file has been made within the input stream. This saves the file for use as a guess of the Hessian matrix in the following calculation. The link can be also done in the shell script.

The generators used to define the C2v symmetry are X and XY, plane yz and axis z. They differ from those used in other examples as in section [*]. The only consequence is that the order of the symmetries in SEWARD differs. In the present case the order is: ${\rm a}_1$, ${\rm a}_2$, ${\rm b}_1$, and ${\rm b}_2$, and consequently the classification by symmetries of the orbitals in the SCF and RASSCF inputs will differ. It is therefore recommended to initially use the option TEST in the SEWARD input to check the symmetry option. This option, however, will stop the calculation after the SEWARD input head is printed.

The calculation converges in three steps. We change now the input. We can choose between replacing by hand the geometry of the SEWARD input or use the same $WorkDir directory and let the program to take the last geometry stored into the communication RUNFILE file. In any case the new input can be:




>>>  Set  Output  File  <<<
>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

>>LINK  $CurrDir/$Project.ForceConstant.STO-3G  COMOLD

  &SEWARD  &END
Title
1,3,-cyclopentadiene  molecule
Symmetry
  X  XY
Basis  set
C.ANO-L...4s3p1d.
C1  .0000000000  .0000000000  -2.3726116671
C2  .0000000000  2.2447443782  -.5623842095
C3  .0000000000  1.4008186026  1.8537195887
End  of  basis
Basis  set
H.ANO-L...2s.
H1  1.6523486260  .0000000000  -3.6022531906
H2  .0000000000  4.1872267035  -1.1903003793
H3  .0000000000  2.5490335048  3.5419847446
End  of  basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<<
  &SCF  &END
TITLE
  cyclopentadiene  molecule
OCCUPIED
9  1  6  2
ITERATIONS
40
END  OF  INPUT

  &RASSCF  &END
LUMORB
TITLE
  cyclopentadiene  molecule  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  9  0  6  0
RAS2
  0  2  0  3
ITER
50,25
CIMX
25
END  OF  INPUT
>>COPY  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<

  &RASSCF  &END
JOBIPH
CIREstart
TITLE
  cyclopentadiene  molecule  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  9  0  6  0
RAS2
  0  2  0  3
ITER
50,25
CIMX
25
END  OF  INPUT
>>COPY  $Project.JobIph  $Project.JobOld

  &ALASKA  &END
End  of  file

  &SLAPAF  &END
OldForce  Constant  Matrix
Iterations
80
Thrs
0.5D-06  1.0D-03
End  of  Input

>>>  EndDo  <<<

The RUNOLD file will be used by SLAPAF as initial Hessian to carry out the relaxation. This use of the RUNFILE can be done between any different calculations provided they work in the same symmetry.

In the new basis set, the resulting optimized geometry at the CASSCF level in C2v symmetry is:



********************************************
* Values of internal coordinates *
********************************************
C2C1 2.851490 Bohr
C3C2 2.545737 Bohr
C3C3 2.790329 Bohr
H1C1 2.064352 Bohr
H2C2 2.031679 Bohr
H3C3 2.032530 Bohr
C1C2C3 109.71 Degrees
C1C2H2 123.72 Degrees
C2C3H3 126.36 Degrees
H1C1H1 107.05 Degrees

Once we have the optimized geometry we can obtain the force field, to compute the force constant matrix and obtain an analysis of the harmonic frequency. This is done by computing the analytical Hessian at the optimized geometry. Notice that this is a single-shot calculation using the MCKINLEY, which will automatically start the MCLR module in case of a frequency calculation.



  &SEWARD  &END
Title
1,3,-cyclopentadiene  molecule
Symmetry
  X  XY
Basis  set
C.ANO-L...4s3p1d.
  C1  0.0000000000  0.0000000000  -2.3483061484
  C2  0.0000000000  2.2245383122  -0.5643712787
  C3  0.0000000000  1.3951643642  1.8424767578
End  of  basis
Basis  set
H.ANO-L...2s.
  H1  1.6599988023  0.0000000000  -3.5754797471
  H2  0.0000000000  4.1615845660  -1.1772096132
  H3  0.0000000000  2.5501642966  3.5149458446
End  of  basis
End  of  Input

  &SCF  &END
TITLE
  cyclopentadiene  molecule
OCCUPIED
9  1  6  2
ITERATIONS
40
END  OF  INPUT
  &RASSCF  &END
TITLE
  cyclopentadiene  molecule  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  9  0  6  0
RAS2
  0  2  0  3
ITER
50,25
CIMX
25
LUMORB
END  OF  INPUT

  &MCKINLEY  &END
Perturbation
  Hessian
End  of  Input

Cyclopentadiene has 11 atoms, that mean 3N = 33 Cartesian degrees of freedom. Therefore the MCLR output will contain 33 frequencies. From those, we are just interested in the 3N-6 = 27 final degrees of freedom that correspond to the normal modes of the system. We will discard from the output the three translational (Ti) and three rotational (Ri) coordinates. The table of characters gives us the classification of these six coordinates: a1 (Tz), a2 (Rz), b2 (Tx,Ry), b1 (Ty,Rx). This information is found in the Seward output:



Character Table for C2v

E s(yz) C2(z) s(xz)
a1 1 1 1 1 z
a2 1 -1 1 -1 xy, Rz, I
b2 1 1 -1 -1 y, yz, Rx
b1 1 -1 -1 1 x, xz, Ry

It is simply to distinguish these frequencies because they must be zero, although and because of numerical un accuracies they will be simply close to zero. In the present calculation the harmonic frequencies, the infrared intensities, and the corresponding normal modes printed below in Cartesian coordinates are the following:



Symmetry a1
==============

1 2 3 4 5 6

Freq. i0.06 847.85 966.07 1044.66 1187.60 1492.41

Intensity: 0.340E-08 0.121E-02 0.533E+01 0.415E+00 0.641E-01 0.393E+01

C1 z 0.12305 0.15517 -0.14426 -0.06780 0.06205 0.02429
C2 y 0.00000 0.19533 0.13649 0.10357 -0.02549 -0.08379
C2 z 0.17402 -0.01781 0.06590 -0.00213 0.03194 -0.06400
C3 y 0.00000 0.02739 -0.06782 0.19744 -0.03160 0.16810
C3 z 0.17402 -0.09921 0.00301 0.07315 -0.09872 -0.03285
H1 x 0.00000 -0.01773 -0.00108 -0.00960 0.00375 0.05160
H1 z 0.17402 0.19434 -0.20604 -0.11326 0.10510 0.12231
H2 y 0.00000 0.15234 0.26633 0.10257 0.14973 0.08294
H2 z 0.17402 -0.15622 0.44925 -0.04930 0.60953 0.48030
H3 y 0.00000 -0.18252 -0.27946 0.49335 0.44520 -0.35864
H3 z 0.17402 0.04882 0.15083 -0.11220 -0.44201 0.34603

7 8 9 10 11

Freq. 1579.76 1633.39 3140.69 3315.45 3341.28

Intensity: 0.473E+01 0.432E+00 0.255E+02 0.143E+02 0.571E+01
...

Symmetry a2
==============

1 2 3 4 5

Freq. i5.81 492.93 663.83 872.55 1235.03

...

Symmetry b2
==============

1 2 3 4 5 6

Freq. i11.15 i0.05 858.71 1020.49 1173.32 1386.18

Intensity: 0.249E-01 0.821E-07 0.259E+01 0.743E+01 0.627E-01 0.163E+00
...
7 8 9 10

Freq. 1424.08 1699.08 3305.25 3334.09

Intensity: 0.966E+00 0.427E+00 0.151E+00 0.302E+02
...

Symmetry b1
==============

1 2 3 4 5 6

Freq. i8.13 0.08 349.36 663.03 881.26 980.60

Intensity: 0.463E-01 0.465E-06 0.505E+01 0.896E+02 0.302E+00 0.169E+02
...

7

Freq. 3159.81

Intensity: 0.149E+02
...

Apart from the six mentioned translational and rotational coordinates There are no imaginary frequencies and therefore the geometry corresponds to a stationary point within the C2v symmetry. The frequencies are expressed in reciprocal centimeters.

After the vibrational analysis the zero-point energy correction and the thermal corrections to the total energy, internal, entropy, and Gibbs free energy. The analysis uses the standard expressions for an ideal gas in the canonical ensemble which can be found in any standard statistical mechanics book. The analysis is performed at different temperatures, for instance:



  Temperature  =  273.00  Kelvin
  ---------------------------

  DeltaG  =  59.90  kcal/mol
  ZPVE  =  60.11  kcal/mol
  TotDeltaU  =  60.72  kcal/mol
  TotDeltaU  -  ZPVE  =  0.61  kcal/mol
  DeltaS  =  3.0001  eu/mol
  U(T&R)  =  1.6275  kcal/mol
  Tot(temp)=  2.2393  kcal/mol

Next, polarizabilities (see below) and isotope shifted frequencies are also displayed in the output.



  ************************************
  *  *
  *  Polarizabilities  *
  *  *
  ************************************



  -34.7624440
  0.0000000  -51.8645093
  0.0000000  0.0000000  -57.7540263

For a graphical representation of the harmonic frequencies one can also use the $Project.freq.molden file as an input to the MOLDEN program.

10.2.2 Excited state optimizations

The calculation of excited states using the ALASKA and SLAPAF codes has no special characteristic. The wave function is defined by the SCF or RASSCF programs. Therefore if we want to optimize an excited state the RASSCF input has to be defined accordingly. It is not, however, an easy task, normally because the excited states have lower symmetry than the ground state and one has to work in low order symmetries if the full optimization is pursued.

Take the example of the thiophene molecule (see fig. [*] in next section). The ground state has C2v symmetry: 1 1A1. The two lowest valence excited states are 21A1 and 11B2. If we optimize the geometries within the C2v symmetry the calculations converge easily for the three states. They are the first, second, and first roots of their symmetry, respectively. But if we want to make a full optimization in C1, or even a restricted one in Cs, all three states belong to the same symmetry representation. The higher the root more difficult is to converge it. A geometry optimization requires single-root optimized CASSCF wave-functions, but, unlike in previous MOLCAS versions, we can now carry out State-Average (SA) CASSCF calculations between different roots. The wave functions we have with this procedure are based on an averaged density matrix, and a further orbital relaxation is required. The MCLR program can perform such a task by means of a perturbational approach. Therefore, if we choose to carry out a SA-CASSCF calculations in the optimization procedure, the Alaska module will automatically start up the MCLR module.

We are going to optimize the three states of thiophene in C2v symmetry. The inputs are:



>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

  &SEWARD  &END
Title
Thiophene  molecule
Symmetry
  X  XY
Basis  set
S.ANO-S...4s3p2d.
S1  .0000000000  .0000000000  -2.1793919255
End  of  basis
Basis  set
C.ANO-S...3s2p1d.
C1  .0000000000  2.3420838459  .1014908659
C2  .0000000000  1.3629012233  2.4874875281
End  of  basis
Basis  set
H.ANO-S...2s.
H1  .0000000000  4.3076765963  -.4350463731
H2  .0000000000  2.5065969281  4.1778544652
End  of  basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<
  &SCF  &END
TITLE
  Thiophene  molecule
OCCUPIED
11  1  7  3
ITERATIONS
40
END  OF  INPUT

  &RASSCF  &END
LUMORB
TITLE
  Thiophene  molecule  1  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  11  0  7  1
RAS2
  0  2  0  3
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<

  &RASSCF  &END
JOBIPH
CIREstart
TITLE
  Thiophene  molecule  1  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  11  0  7  1
RAS2
  0  2  0  3
ITER
50,25
END  OF  INPUT

  &ALASKA  &END
End  of  Input

  &SLAPAF  &END
Iterations
20
Thrs
0.5D-06  1.0D-03
End  of  Input
>>>  ENDDO  <<<

for the ground state. For the two excited states we will replace the RASSCF inputs with



  &RASSCF  &END
LUMORB
*JOBIPH
*CIRESTART
TITLE
  Thiophene  molecule  2  1A1
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  11  0  7  1
RAS2
  0  2  0  3
CIROOT
2  2
1  2
1  1
LEVSHFT
1.0
ITER
50,25
RLXRoot
2
END  OF  INPUT

for the 21A1 state. Notice that we are doing a SA-CASSCF calculation including two roots, therefore we must use the keyword RLXROOT within the RASSCF input to specify for which state we want the root. We have also



  &RASSCF  &END
LUMORB
*JOBIPH
*CIRESTART
TITLE
  Thiophene  molecule  1  1B2
SYMMETRY
  2
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  11  0  7  1
RAS2
  0  2  0  3
LEVSHFT
1.0
ITER
50,25
END  OF  INPUT

for the 11B2 state.

To help the program to converge we can include one or more initial RASSCF inputs in the input file. The following is an example for the calculation of the of the 31A' state of thiophene (Cs symmetry) with a previous calculation of the ground state to have better starting orbitals. The option SALA equal to three is used to relax the CASSCF orbitals for the exact root which we are interested.



>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

  &SEWARD  &END
Title
Thiophene  molecule
Symmetry
  X
Basis  set
S.ANO-S...4s3p2d.
S1  .0000000000  .0000000000  -2.1174458547
End  of  basis
Basis  set
C.ANO-S...3s2p1d.
C1  .0000000000  2.4102089951  .1119410701
C1b  .0000000000  -2.4102089951  .1119410701
C2  .0000000000  1.3751924147  2.7088559532
C2b  .0000000000  -1.3751924147  2.7088559532
End  of  basis
Basis  set
H.ANO-S...2s.
H1  .0000000000  4.3643321746  -.4429940876
H1b  .0000000000  -4.3643321746  -.4429940876
H2  .0000000000  2.5331491787  4.3818833166
H2b  .0000000000  -2.5331491787  4.3818833166
End  of  basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<
  &SCF  &END
TITLE
  Thiophene  molecule
OCCUPIED
18  4
ITERATIONS
40
END  OF  INPUT

  &RASSCF  &END
LUMORB
TITLE
  Thiophene  molecule  1A'
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  18  1
RAS2
  0  5
CIROOT
1  1
1
LEVSHFT
1.0
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld

  &RASSCF  &END
JOBIPH
TITLE
  Thiophene  molecule  3  1A'
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  18  1
RAS2
  0  5
CIROOT
3  3
1  2  3
1  1  1
LEVSHFT
1.0
ITER
50,25
RLXRoot
3
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<


  &RASSCF  &END
JOBIPH
CIRESTART
TITLE
  Thiophene  molecule  3  1A'
SYMMETRY
  1
SPIN
  1
NACTEL
  6  0  0
INACTIVE
  18  1
RAS2
  0  5
CIROOT
3  3
1  2  3
1  1  1
LEVSHFT
1.0
ITER
50,25
RLXRoot
3
END  OF  INPUT

  &ALASKA  &END
End  of  Input

  &SLAPAF  &END
End  of  Input
>>>  ENDDO  <<<

It should be remembered that geometry optimizations for excited states are difficult. Not only can it be difficult to converge the corresponding RASSCF calculation, but we must also be sure that the order of the states does not change during the optimization of the geometry. This is not uncommon and the optimization must be followed by the user.

Sometimes may be interesting to follow the path of the optimization by looking at each one of the output files generated by MOLCASAll the iterative information is stored in the input file if the "Set Output File" command as not used. If it was used the output files of each complete iteration are stored in the $WorkDir directory under the names 1.save.$iter, for instance: 1.save.1, 1.save.2, etc. You should not remove the $WorkDir directory if you want to keep them.

10.2.3 Restrictions in symmetry or geometry.


10.2.3.1 Optimizing with geometrical constraints.

A common situation in geometry optimizations is to have one or several coordinates fixed or constrained and vary the remaining coordinates. As an example we will take the biphenyl molecule, two benzene moieties bridged by a single bond. The ground state of the molecule is not planar. One benzene group is twisted by 44$\,^o$ degrees with respect to the other [210]. We can use this example to perform two types of restricted optimizations. The simplest way to introduce constraints is to give a coordinate a fixed value and let the other coordinates to be optimized. For instance, let's fix the dihedral angle between both benzenes to be fixed to 44$\,^o$ degrees. Within this restriction, the remaining coordinates will be fully optimized. The Constraints keyword in the program SLAPAF will take care of the restriction. The input could be:



>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<
  &SEWARD  &END
Title
  Biphenyl  twisted  D2
Symmetry
  XY  XZ
Basis  set
C.ANO-S...3s2p1d.
C1  1.4097582886  .0000000000  .0000000000
C2  2.7703009377  2.1131321616  .8552434921
C3  5.4130377085  2.1172148045  .8532344474
C4  6.7468359904  .0000000000  .0000000000
End  of  basis
Basis  set
H.ANO-S...2s.
H2  1.7692261798  3.7578798540  1.5134152112
H3  6.4188773347  3.7589592975  1.5142479153
H4  8.7821560635  .0000000000  .0000000000
End  of  basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<
  &SCF  &END
TITLE
  Biphenyl  twisted  D2
OCCUPIED
  12  9  9  11
ITERATIONS
50
END  OF  INPUT

  &RASSCF  &END
LUMORB
TITLE
  Biphenyl  twisted  D2
SYMMETRY
  1
SPIN
  1
NACTEL
  12  0  0
INACTIVE
  11  7  7  10
RAS2
  2  4  4  2
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<

  &RASSCF  &END
JOBIPH
CIRESTART
TITLE
  Biphenyl  twisted  D2
SYMMETRY
  1
SPIN
  1
NACTEL
  12  0  0
INACTIVE
  11  7  7  10
RAS2
  2  4  4  2
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld

  &ALASKA  &END
End  of  input

  &SLAPAF  &END
Constraints
d1  =  Dihedral  C2  C1  C1(XY)  C2(XY)
Values
d1  =  44.4  degrees
End  of  Constraints
Iterations
30
End  of  Input
>>>  ENDDO  <<<

One important consideration about the constraint. You do not need to start at a geometry having the exact value for the coordinate you have selected (44.4 degrees for the dihedral angle here). The optimization will lead you to the right solution. On the other hand, if you start exactly with the dihedral being 44.4 deg the code does not necessarily will freeze this value in the first iterations, but will converge to it at the end. Therefore, it may happen that the value for the dihedral differs from the selected value in the initial iterations. You can follow the optimization steps in the $WorkDir directory using the MOLDEN files generated automatically by MOLCAS.

Now we will perform the opposite optimization: we want to optimize the dihedral angle relating both benzene units but keep all the other coordinates fixed. We could well use the same procedure as before adding constraints for all the remaining coordinates different from the interesting dihedral angle, but to build the input would be tedious. Therefore, instead of keyword Constraints we will make use of the keywords Vary and Fix.

The input file should be:



>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

  &SEWARD  &END
Title
  Biphenyl  twisted  D2
Symmetry
  XY  XZ
Basis  set
C.ANO-S...3s2p1d.
C1  1.4097582886  .0000000000  .0000000000
C2  2.7703009377  2.1131321616  .8552434921
C3  5.4130377085  2.1172148045  .8532344474
C4  6.7468359904  .0000000000  .0000000000
End  of  basis
Basis  set
H.ANO-S...2s.
H2  1.7692261798  3.7578798540  1.5134152112
H3  6.4188773347  3.7589592975  1.5142479153
H4  8.7821560635  .0000000000  .0000000000
End  of  basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<
  &SCF  &END
TITLE
  Biphenyl  twisted  D2
OCCUPIED
  12  9  9  11
ITERATIONS
50
END  OF  INPUT

  &RASSCF  &END
LUMORB
TITLE
  Biphenyl  twisted  D2
SYMMETRY
  1
SPIN
  1
NACTEL
  12  0  0
INACTIVE
  11  7  7  10
RAS2
  2  4  4  2
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<

  &RASSCF  &END
JOBIPH
CIRESTART
TITLE
  Biphenyl  twisted  D2
SYMMETRY
  1
SPIN
  1
NACTEL
  12  0  0
INACTIVE
  11  7  7  10
RAS2
  2  4  4  2
ITER
50,25
END  OF  INPUT
!cp  $Project.JobIph  $Project.JobOld

  &ALASKA  &END
End  of  input

  &SLAPAF  &END
Internal  coordinates
b1  =  Bond  C1  C1(XY)
b2  =  Bond  C1  C2
b3  =  Bond  C2  C3
b4  =  Bond  C3  C4
h1  =  Bond  C2  H2
h2  =  Bond  C3  H3
h3  =  Bond  C4  H4
a1  =  Angle  C2  C1  C1(XY)
a2  =  Angle  C1  C2  C3
a3  =  Angle  C1  C2  H2
a4  =  Angle  C2  C3  H3
phi  =  Dihedral  C2  C1  C1(XY)  C2(XY)
d1  =  Dihedral  H2  C2  C1  C1(XY)
d2  =  OutOfP  C3  C1(XY)  C1  C2
d3  =  Dihedral  H3  C3  C2  H2
Vary
  phi  =  1.0  phi
Fix
  b1  =  1.0  b1
  b2  =  1.0  b2
  b3  =  1.0  b3
  b4  =  1.0  b4
  h1  =  1.0  h1
  h2  =  1.0  h2
  h3  =  1.0  h3
  a1  =  1.0  a1
  a2  =  1.0  a2
  a3  =  1.0  a3
  a4  =  1.0  a4
  d1  =  1.0  d1
  d2  =  1.0  d2
  d3  =  1.0  d3
End  of  Internal
Iterations
30
End  of  Input
>>>  ENDDO  <<<

To be able to optimize the molecule in that way a D2 symmetry has to be used. In the definition of the internal coordinates we can use an out-of-plane coordinate: C2 C2(xy) C1(xy) C1 or a dihedral angle C2 C1 C1(xy) C2(xy). In this case there is no major problem but in general one has to avoid as much as possible to define dihedral angles close to 180$\,^o$ ( trans conformation ). The SLAPAF program will warn about this problem if necessary. In the present example, angle 'phi' is the angle to vary while the remaining coordinates are frozen. All this is only a problem in the user-defined internal approach, not in the non-redundant internal approach used by default in the program. In case we do not have the coordinates from a previous calculation we can always run a simple calculation with one iteration in the SLAPAF program.

It is not unusual to have problems in the relaxation step when one defines internal coordinates. Once the program has found that the definition is consistent with the molecule and the symmetry, it can happen that the selected coordinates are not the best choice to carry out the optimization, that the variation of some of the coordinates is too large or maybe some of the angles are close to their limiting values ($\pm$180$\,^o$ for Dihedral angles and $\pm$90$\,^o$ for Out of Plane angles). The SLAPAF program will inform about these problems. Most of the situations are solved by re-defining the coordinates, changing the basis set or the geometry if possible, or even freezing some of the coordinates. One easy solution is to froze this particular coordinate and optimize, at least partially, the other as an initial step to a full optimization. It can be recommended to change the definition of the coordinates from internal to Cartesian.

Figure 10.3: Twisted biphenyl molecule
\begin{figure}{---------------------------------------------------}
\myincludegraphics{advanced.examples/biphenyl}\end{figure}


10.2.3.2 Optimizing with symmetry restrictions.

Presently, MOLCAS is prepared to work in the point groups C1, Ci, Cs, C2, D2, C2h, C2v, and D2h. To have the wave functions or geometries in other symmetries we have to restrict orbital rotations or geometry relaxations specifically. We have shown how to in the RASSCF program by using the SUPSym option. In a geometry optimization we may also want to restrict the geometry of the molecule to other symmetries. For instance, to optimize the benzene molecule which belongs to the D6h point group we have to generate the integrals and wave function in D2h symmetry, the highest group available, and then make the appropriate combinations of the coordinates chosen for the relaxation in the SLAPAF program, as is shown in the manual.

As an example we will take the ammonia molecule, NH3. There is a planar transition state along the isomerization barrier between two pyramidal structures. We want to optimize the planar structure restricted to the D3h point group. However, the electronic wave function will be computed in Cs symmetry (C2v is also possible) and will not be restricted, although it is possible to do that in the RASSCF program.

The input for such a geometry optimization is:



>>>  Set  MaxIter  50  <<<
>>>  Do  while  <<<

  &SEWARD  &END
Title
  NH3,  planar
Symmetry
  Z
Basis  Set
N.ANO-L...4s3p2d.
N  .0000000000  .0000000000  .0000000000
End  of  Basis
Basis  set
H.ANO-L...3s2p.
H1  1.9520879910  .0000000000  .0000000000
H2  -.9760439955  1.6905577906  .0000000000
H3  -.9760439955  -1.6905577906  .0000000000
End  of  Basis
End  of  Input

>>>  IF  (  ITER  =  1  )  <<<
  &SCF  &END
Title
  NH3,  planar
Occupied
  4  1
Iterations
40
End  of  Input

  &RASSCF  &END
LUMORB
Title
  NH3,  planar
Symmetry
1
Spin
1
Nactel
  8  0  0
INACTIVE  ORBITALS
  1  0
RAS2  ORBITALS
  6  2
ITER
50,20
End  of  Input
!cp  $Project.JobIph  $Project.JobOld
>>>  ENDIF  <<<

  &RASSCF  &END
JOBIPH
CIRESTART
Title
  NH3,  planar
Symmetry
1
Spin
1
Nactel
  8  0  0
INACTIVE  ORBITALS
  1  0
RAS2  ORBITALS
  6  2
ITER
50,20
End  of  Input
!cp  $Project.JobIph  $Project.JobOld

  &ALASKA  &END
End  of  input

  &SLAPAF  &END
Internal  coordinates
b1  =  Bond  N  H1
b2  =  Bond  N  H2
b3  =  Bond  N  H3
a1  =  Angle  H1  N  H2
a2  =  Angle  H1  N  H3
Vary
r1  =  1.0  b1  +  1.0  b2  +  1.0  b3
Fix
r2  =  1.0  b1  -  1.0  b2
r3  =  1.0  b1  -  1.0  b3
a1  =  1.0  a1
a2  =  1.0  a2
End  of  internal
Iterations
20
End  of  input
>>>  ENDDO  <<<

All four atoms are in the same plane. Working in Cs, planar ammonia has five degrees of freedom. Therefore we must define five independent internal coordinates, in this case the three N-H bonds and two of the three angles H-N-H. The other is already defined knowing the two other angles. Now we must define the varying coordinates. The bond lengths will be optimized, but all three N-H distances must be equal. First we define (see definition in the previous input) coordinate r1 equal to the sum of all three bonds; then, we define coordinates r2 and r3 and keep them fixed. r2 will ensure that bond1 is equal to bond2 and r3 will assure that bond3 is equal to bond1. r2 and r3 will have a zero value. In this way all three bonds will have the same length. As we want the system constrained into the D3h point group, the three angles must be equal with a value of 120 degrees. This is their initial value, therefore we simply keep coordinates ang1 and ang2 fixed. The result is a D3h structure:



********************************************
* Values of primitive internal coordinates *
********************************************
~
B1 : Bond Length= 1.8957/ bohr
B2 : Bond Length= 1.8957/ bohr
B3 : Bond Length= 1.8957/ bohr
A1 : Angle= 120.0000/degree, 2.0944/rad
A2 : Angle= 120.0000/degree, 2.0944/rad

In a simple case like this an optimization without restrictions would also end up in the same symmetry as the initial input.


10.2.4 Optimizing with Z-Matrix.

An alternative way to optimize a structure with geometrical and/or symmetrical constraints is to combine the Z-Matrix definition of the molecular structure used for the program SEWARD with a coherent definition for the Internal Coordinated used in the optimization by program SLAPAF.

Here is an examples of optimization of the methyl carbanion. Note that the wavefunction is calculated within the Cs symmetry but the geometry is optimized within the C3v symmetry through the ZMAT and the Internal Coordinates definitions.



>>> Set MaxIter 500 <<<
>>> Do While <<<

&SEWARD &END
Symmetry
Y

ZMAT
H.Aug-cc-pVDZ....
C.Aug-cc-pVDZ....

C1
X2 1 1.00
H3 1 1.09 2 105.
H4 1 1.09 2 105. 3 120.

End of input

&SCF &END
Charge
-1
End of input

&ALASKA &END
End of input

&SLAPAF &END
Internal Coordinates
CX2 = Bond C1 X2
CH3 = Bond C1 H3
CH4 = Bond C1 H4
XCH3 = Angle X2 C1 H3
XCH4 = Angle X2 C1 H4
DH4 = Dihedral H3 X2 C1 H4
Vary
SumCH34 = 1. CH3 +2. CH4
SumXCH34 = 1. XCH3 +2. XCH4
Fix
rCX2 = 1.0 CX2
DifCH34 = 2. CH3 -1. CH4
DifXCH34 = 2. XCH3 -1. XCH4
dDH4 = 1.0 DH4
End of Internal
PRFC
Iterations
10
End of input

>>> EndDo <<<

Note that the dummy atom X2 is used to define the Z axis and the planar angles for the hydrogen atoms. The linear combinations of bond distances and planar angles in the expression in the Vary and Fix sections are used to impose the C3v symmetry.

Another examples where the wavefunction and the geometry can be calculated within different symmetry groups is benzene. In this case, the former uses D2h symmetry and the latter D6h symmetry. Two special atoms are used: the dummy X1 atom defines the center of the molecule while the ghost Z2 atom is used to define the C6 rotational axis (and the Z axis).



>>> Set MaxIter 500 <<<
>>> Do While <<<

&SEWARD &END
Symmetry
X Y Z

ZMAT
H.ANO-S...2s.
C.ANO-S...3s2p.

X1
Z2 1 1.00
C3 1 1.3915 2 90.
C4 1 1.3915 2 90. 3 60.
H5 1 2.4715 2 90. 3 0.
H6 1 2.4715 2 90. 3 60.

End of input

&SCF &END
End of input

&ALASKA &END
End of input

&SLAPAF &END
Internal Coordinates
XC3 = Bond X1 C3
XC4 = Bond X1 C4
XH5 = Bond X1 H5
XH6 = Bond X1 H6
CXC = Angle C3 X1 C4
HXH = Angle H5 X1 H6
Vary
SumC = 1.0 XC3 + 2.0 XC4
SumH = 1.0 XH5 + 2.0 XH6
Fix
DifC = 2.0 XC3 - 1.0 XC4
DifH = 2.0 XH5 - 1.0 XH6
aCXC = 1.0 CXC
aHXH = 1.0 HXH
End of Internal
PRFC
Iterations
10
End of input

>>> EndDo <<<

Note that the ghost atom Z2 is used to define the geometry within the Z-Matrix but it does not appear in the Internal Coordinates section. On the other hand, the dummy atom X1 represents the center of the molecule and it is used in the Internal Coordinates section.


10.2.5 CASPT2 optimizations

For systems showing a clear multiconfigurational nature, the CASSCF treatment on top of the HF results is of crucial importance in order to recover the large non dynamical correlation effects. On the other hand, ground-state geometry optimizations of closed shell systems are not exempt from non dynamical correlation effects. In general, molecules containing $\pi$-electrons suffer from significant effects of non dynamical correlation, even more in presence of conjugated groups. Several studies on systems with delocalized bonds have shown the effectiveness of the CASSCF approach in reproducing the main geometrical parameters with high accuracy [211,212,213].

However, pronounced effects of dynamical correlation often occur in systems with $\pi$-electrons, especially in combination with polarized bonds. An example is given by the C=O bond length, which is known to be very sensitive to an accurate description of the dynamical correlation effects [214]. We will show now that the inherent limitations of the CASSCF method can be successfully overcome by employing a CASPT2 geometry optimization, which uses a numerical gradient procedure of recent implementation. A suitable molecule for this investigation is acrolein. As many other conjugated aldehydes and ketones, offers an example of s-cis/s-trans isomerism (Figure [*]). Due to the resonance between various structures involving $\pi$ electrons, the bond order for the C-C bond is higher than the one for a non-conjugated C-C single bond. This partial double-bond character restricts the rotation about such a bond, giving rise to the possibility of geometrical isomerism, analogue to the cis-trans one observed for conventional double bonds.

A CASPT2 geometry optimization in MOLCAS can be performed using AUTO program. A possible input for the CASPT2 geometry optimization of the s-trans isomer is displayed below. The basis set is ridiculously small, making the use of CASPT2 almost meaningless: Change to, e.g., ANO-RCC-VDZP of larger for a serious run. But remember that small, quick calculations have more uses than merely checking that the job script works. For example, the small calculation can be done using just RASSCF, the runfile can be saved for use in a subsequent structure optimization with large basis at the CASPT2 level, with starting geometry already close to the correct ones. Also, in other cases than this, it is valuable to have the small-basis orbital file expanded (using EXPBAS) to provide proper starting orbitals for the RASSCF. Starting with SCF orbitals may have orbitals in an improper order for getting the active orbitals right.



>>>  SET  MAXITER  500  <<<

&GATEWAY
Title
  Acrolein  Cs  symmetry  -  transoid
Coord
  8

  O  0.0000000  -1.4511781  -1.3744831
  C  0.0000000  -0.8224882  -0.1546649
  C  0.0000000  0.7589531  -0.0387200
  C  0.0000000  1.3465057  1.2841925
  H  0.0000000  -1.4247539  0.8878671
  H  0.0000000  1.3958142  -1.0393956
  H  0.0000000  0.6274048  2.2298215
  H  0.0000000  2.5287634  1.4123985
Group
  X
*  Demonstration  run  with  very  small  basis,  making  CASPT2  almost  meaningless.
*  Change  to  e.g.  ANO-RCC-VDZP  for  a  more  serious  run.
Basis
ANO-S-MB
End  of  Input

>>>>>>>>>>>>>  Do  while  <<<<<<<<<<<<

&SEWARD
End  of  Input

>>>>>>>>  IF  (  ITER  =  1  )  <<<<<<<<<<<
&SCF
Title
Acrolein  Cs  symmetry
Occupied
*The  symmetry  species  are  a'  a''
  13  2
End  of  input
>>>>>>>  ENDIF  <<<<<<<<<<<<<<<<<<<<<

&RASSCF
LUMORB
Title
  Acrolein  ground  state
nActEl
  4  0  0
Inactive
  13  0
*The  symmetry  species  are  a'  a''
Ras2
  0  4
End  of  input

&CASPT2
End  of  input

&SLAPAF
End  of  input

>>>>>>>>>>>>>  ENDDO  <<<<<<<<<<<<<<

Experimental investigations assign a planar structure for both the isomers. We can take advantage of this result and use a Cs symmetry throughout the optimization procedure. Moreover, the choice of the active space is suggested by previous calculations on analogous systems. The active space contains 4 $\pi$ MOs /4 $\pi$ electrons, thus what we will call shortly a $\pi$-CASPT2 optimization.

The structure of the input follows the trends already explained in other geometry optimizations, that is, loops over the set of programs ending with SLAPAF. Notice that CASPT2 optimizations require obviously the CASPT2 input, but also the input for the ALASKA program, which computes the gradient numerically. Apart from that, a CASPT2 optimization input is identical to the corresponding CASSCF input. With the present implementation of the CASPT2 numerical gradients procedure, the output of our calculation will contain by default all the intermediate outputs at each point of the grid. Even for a calculation with a small number of internal coordinates, slowly convergent situations can end up in an undesirably huge output file. To avoid this annoying circumstance, the keyword SET OUTPUT FILE of the AUTO script can be used. In this case, the output will be split into several files each corresponding the the whole calculations done in a single iteration of the geometry optimization. The keyword SET OVER can be also useful especially when we want to save disk space. It replaces the output of each iteration with the one of the subsequent, overriding the previous output file. The output file corresponding to the final iteration, contains the summary concerning each iterations. For our example the résumé is the following:




************************************************************************************************
***************** Energy Statistics for Geometry Optimization **********************************
************************************************************************************************

Energy Grad Grad Step Estimated Hess Geom Hess
Iter Energy Change Norm Max Element Max Element Final Ene Index Upd Update
1 -191.46171057 0.00000000 0.028786 0.015590 nrc001 0.019477 nrc001 -191.46189776 0 RF(S) None
2 -191.46195237-0.00024180 0.010010 0.005154 nrc001 0.012652 nrc001 -191.46198941 0 RF(S) BFGS
3 -191.46198838-0.00003601 0.001320-0.000670 nrc012 -0.002159 nrc012 -191.46198994 0 RF(S) BFGS
4 -191.46199164-0.00000326 0.000661-0.000265 nrc012 0.000527 nrc004 -191.46199187 0 RF(S) BFGS
5 -191.46199182-0.00000017 0.000435-0.000184 nrc012 -0.000903 nrc012 -191.46199231 0 RF(S) BFGS
6 -191.46199215-0.00000034 0.000190 0.000094 nrc010 0.000982 nrc010 -191.46199230 0 RF(S) BFGS
7 -191.46199217-0.00000002 0.000086 0.000028 nrc010 0.000570 nrc012 -191.46199220 0 RF(S) BFGS
8 -191.46199214 0.00000003 0.000134 0.000061 nrc010 0.000464 nrc010 -191.46199217 0 RF(S) BFGS
9 -191.46199216-0.00000002 0.000079 0.000028 nrc010 -0.000232 nrc010 -191.46199215 0 RF(S) BFGS


Cartesian Displacements Gradient in internals
Value Threshold Converged? Value Threshold Converged?
+----------------------------------+----------------------------------+
RMS + 0.2558E-03 0.4000E-03 Yes + 0.2291E-04 0.1000E-03 Yes +
+----------------------------------+----------------------------------+
Max + 0.3879E-03 0.6000E-03 Yes + 0.2761E-04 0.1500E-03 Yes +
+----------------------------------+----------------------------------+

Geometry is converged

************************************************************************************************
************************************************************************************************

The calculation converges in 9 iterations. At this point it is worth noticing how the convergence of CASPT2 energy is not chosen among the criteria for the convergence of the structure. The final structure is in fact decided by checking the Cartesian displacements and the gradient in non-redundant internal coordinates.

CASPT2 optimizations are expensive. Notice that they are based on numerical gradients and many point-wise calculations are needed. In particular, double-sided gradients are computed in Cartesian. Therefore, each macro-iteration in the optimization requires 2*N + 1 Seward/RASSCF/CASPT2 calculations, with N being the Cartesian degrees of freedom. In the present example, acrolein has eight atoms. From each atom, only two Cartesian coordinates are free to move (we are working within the Cs symmetry and the third coordinate is frozen), therefore the total number of Seward/RASSCF/CASPT2 iterations within each macro-iteration is 2*(8*2) + 1, that is, 33. It is not an easy task.

The Table [*] displays the equilibrium geometrical parameters computed at the $\pi$-CASSCF and $\pi$-CASPT2 level of theory for the ground state of both isomers of acrolein. For sake of comparison, Table [*] includes experimental data obtained from microwave spectroscopy studies[215]. The computed parameters at $\pi$-CASPT2 level are in remarkable agreement with the experimental data. The predicted value of the C=C bond length is very close to the double bond length observed in ethylene. The other C-C bond has a length within the range expected for a C-C single bond: it appears shorter in the s-trans isomer as a consequence of the reduction of steric hindrance between the ethylenic and aldehydic moieties. CASSCF estimates a carbon-oxygen bond length shorter than the experimental value. For $\pi$-CASSCF optimization in conjugated systems this can be assumed as a general behavior [216,214]. To explain such a discrepancy, one may invoke the fact that the C=O bond distance is particularly sensitive to electron correlation effects. The $\pi$ electron correlation effects included at the $\pi$-CASSCF level tend to overestimate bond lengths. However, the lack of $\sigma$ electron correlation, goes in the opposite direction, allowing shorter bond distances for double bonds. For the C-C double bonds, these contrasting behaviors compensate each other [213] resulting in quite an accurate value for the bond length at the $\pi$-CASSCF level. On the contrary, the extreme sensitivity of the C=O bond length to the electron correlation effects, leads to a general underestimation of the C-O double bond lengths, especially when such a bond is part of a conjugated system. It is indeed the effectiveness of the CASPT2 method in recovering dynamical correlation which leads to a substantial improvement in predicting the C-O double bond length.

\begin{figure}{---------------------------------------------------}
\end{figure}
\myincludegraphics{advanced.examples/acrolein}


Table 10.8: Geometrical parameters for the ground state of acrolein
Parametersa $\pi$-CASSCF [04/4]   $\pi$-CASPT2   Expt.b
  s-cis s-trans   s-cis s-trans    
C1=O 1.204 1.204   1.222 1.222   1.219
C1-C2 1.483 1.474   1.478 1.467   1.470
C2=C3 1.340 1.340   1.344 1.344   1.345
$\angle$ C1C2C3 123.0 121.7   121.9 120.5   119.8
$\angle$ C2C1O 124.4 123.5   124.5 124.2   -
aBond distances in Å and angles in degrees.
bMicrowave spectroscopy data from ref. [215]. No difference between s-cis and s-trans isomers is reported

The use of numerical CASPT2 gradients can be extended to all the optimizations available in SLAPAF, for instance transition state searches. Use the following input for the water molecule to locate the linear transition state:



&GATEWAY
Title
  Water,  STO-3G  Basis  set
Coord
3

H1  -0.761622  0.000000  -0.594478
H2  0.761622  0.000000  -0.594478
O  0.000000  0.000000  0.074915
Basis  set
  STO-3G
Group
  NoSym
End  of  input

>>>  SET  MAXITER  500  <<<
>>  DO  WHILE

&SEWARD
End  of  input

>>>  IF  (  ITER  =  1  )  <<<
&SCF
Title
water,  STO-3g  Basis  set
Occupied
5
End  of  input
>>>  ENDIF  <<<

&RASSCF
LumOrb
Nactel
2  0  0
Inactive
4
Ras2
2
End  of  Input

&CASPT2
End  of  input

&SLAPAF
TS
End  of  Input
>>>  ENDDO  <<<

After ten macro-iterations the linear water is reached:



**********************************************************************************************
*************** Energy Statistics for Geometry Optimization **********************************
**********************************************************************************************

Energy Grad Grad Step Estimated Hess Geom Hess
Iter Energy Change Norm Max Element Max Element Final Ener Index Upd Upd
1 -75.00603925 0.00000000 0.000505-0.000333 nrc001 0.149153 nrc003 -75.00437273 1 MFRFS None
2 -75.00256314 0.00347612 0.033114-0.027092 nrc003 0.145257 nrc003 -75.00034225 1 MFRFS MSP
3 -74.99310559 0.00945755 0.083775-0.078714 nrc003 -0.184679 nrc002 -74.98597057 1 MFRFS MSP
4 -74.97219951 0.02090608 0.163015 0.086748 nrc002 0.226701 nrc003 -74.94503565 1 MFRFS MSP
5 -74.93277784 0.03942168 0.201340 0.123100 nrc002 0.223518 nrc003 -74.88089265 1 MFRFS MSP
6 -74.89601350 0.03676433 0.148875-0.100716 nrc003 0.230650 nrc003 -74.86658241 1 MFRFS MSP
7 -74.87796405 0.01804946 0.044034 0.037502 nrc002 0.055719 nrc002 -74.87855019 1 MFRFS MSP
8 -74.87878116-0.00081712 0.009948-0.007364 nrc003 0.016283 nrc003 -74.87873039 1 MFRFS MSP
9 -74.87872319 0.00005797 0.000848-0.000641 nrc001 0.001232 nrc003 -74.87872320 1 MFRFS MSP
10 -74.87872373-0.00000054 0.000053-0.000048 nrc001 -0.000066 nrc001 -74.87872373 1 MFRFS MSP


Cartesian Displacements Gradient in internals
Value Threshold Converged? Value Threshold Converged?
+----------------------------------+----------------------------------+
RMS + 0.4174E-04 0.1200E-02 Yes + 0.3719E-04 0.3000E-03 Yes +
+----------------------------------+----------------------------------+
Max + 0.4665E-04 0.1800E-02 Yes + 0.4836E-04 0.4500E-03 Yes +
+----------------------------------+----------------------------------+

Geometry is converged

***********************************************************************************************
***********************************************************************************************


next up previous contents
Next: 10.3 Computing a reaction Up: 10. Examples Previous: 10.1 Computing high symmetry