To optimize a molecular geometry is probably one of the most frequent interests of a quantum chemist [16]. In the present section we examine some examples of obtaining stationary points on the energy surfaces by using the analytical derivative program ALASKA which computes the gradient for optimized wave functions. This is presently available at the SCF and RASSCF levels of calculation. 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. MOLCAS-4 simplifies the task with respect to previous releases because there is no need 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 coodinates and displacements [17]. In any case, in the following examples we are going to keep the definition of the internal coordinates. As they depend on the symmetry of the system it might be somewhat difficult in some systems to define them. It is, however, strongly recommended to let the program define its own set of non-redundant internal coordinates.
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 numerical hessian calculations to compute the force fields at the optimized geometries.
We begin by defining the input for the initial calculation.
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 
symmetry.
We will use the following input where the basis set has been introduced directly in the SEWARD input:
&STRUCTURE &END
Method casscf
!ln -fs $TempDir/$Project.Relax.STO-3G RELAX
End of input
&SEWARD &END
Title
1,3,-cyclopentadiene. STO-3G basis set.
Symmetry
X XY
Basis set
C....2s1p. / Inline
6.0 1
6 2
71.616837
13.045096
3.530512
2.941249
0.683483
0.222290
0.154329 0.0
0.535328 0.0
0.444635 0.0
0.0 -0.099967
0.0 0.399513
0.0 0.700115
3 1
2.941249
0.683483
0.222290
0.155916
0.607684
0.391957
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....1s. / Inline
1.0 0
3 1
3.425251
0.623914
0.168855
0.154329
0.535328
0.444635
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
NoInvariance
End of Input
&SLAPAF &END
Internal coordinates
b1 = Bond C1 C2
b2 = Bond C2 C3
b3 = Bond C3 C3(xy)
b4 = Bond C1 H1
b5 = Bond C2 H2
b6 = Bond C3 H3
a1 = Angle C1 C2 C3
a2 = Angle C1 C2 H2
a3 = Angle C2 C3 H3
a4 = Angle H1 C1 H1(x)
Vary
b1
b2
b3
b4
b5
b6
a1
a2
a3
a4
End of Internal
Iterations
40
Thrs
0.5D-06 1.0D-03
End of Input
A link to the RELAX file has been made within the STRUCTURE namelist. 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
symmetry are X and XY, plane 
and axis 
. They
differ from those used in other examples as in section 4.1.1.
The only consequence is that the order of the symmetries in SEWARD
differs. In the present case the order is: 
, 
, 
, and 
,
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.
Within symmetry some of the internal coordinates of the
whole molecule are redundant. We have to specify only the distinct
internal coordinates. The molecule is placed in the 
plane
along the axis as the binary axis (see Figure 4.1).
The carbon atoms can be defined as C1, C2, C2(xy), C3, and C3(xy), where we have
in the parenthesis the symmetry operation needed to relate the
atom with its equivalent center. C1 is unique and is placed in
the axis. From a total of five
carbon-carbon bond lengths in the molecule we must define
three. Bond C1-C2, bond C2-C3, and bond C3-C3(xy). The other
two bonds C1-C2(xy) and C2(xy)-C3(xy) are already defined by
the symmetry of the molecule.
As regards the angles, the ring has three distinct angles
among the carbons, but, once we have defined the bonds,
only the angle C1-C2-C3 is necessary. Angle C2-C1-C2(xy) is
defined once the bond distances C1-C2, C2-C3 and the
angle C1-C2-C3 are defined, and angle C2-C3-C3(xy) is already
defined by the angle C1-C2-C3 and the distance C3-C3(xy).
The position of H1 and its equivalent atom with negative 
value can be defined in different ways but all of them
must keep the symmetry. One possibility is
Angle H1 C1 H1(x), but it is also possible to use
Dihedral C3 C2 C1 H1. In both cases H1 and H1(x) are restricted
to axis .
Once the calculation has converged, we replace the STRUCTURE input with:
&STRUCTURE &END Method casscf !ln -fs $TempDir/$Project.Relax.STO-3G RLXOLD End of input
and the SEWARD input by increasing the basis set and using the optimized geometry:
&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
and add to the SLAPAF input the keyword:
OldForce Constant Matrix
The RLXOLD file will be used by SLAPAF as initial Hessian to carry out the relaxation. This use of the RELAX 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 symmetry is:
******************************************** * Values of primitive internal coordinates * ******************************************** B1 : Bond Length= 2.8513/ bohr B2 : Bond Length= 2.5463/ bohr B3 : Bond Length= 2.7908/ bohr B4 : Bond Length= 2.0644/ bohr B5 : Bond Length= 2.0319/ bohr B6 : Bond Length= 2.0328/ bohr A1 : Angle= 109.7203/degree, 1.9150/rad A2 : Angle= 123.8103/degree, 2.1609/rad A3 : Angle= 126.3077/degree, 2.2045/rad A4 : Angle= 106.8958/degree, 1.8657/rad
using the previous internal coordinates definition.
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 numerical Hessian at the optimized geometry.
If we continue in
we will only obtain the analysis of the totally symmetric
modes ( in this case). If we want a full analysis we must
remove all symmetry, that is, work in C
symmetry.
Computing a numerical Hessian is an expensive task because
it takes 2
+1 iterations, where is the number of degrees
of freedom of the molecule within the symmetry. For cyclopentadiene
we have a total of 27 degrees of freedom in C, what means 55 iterations.
To compute simply the in-plane vibrational modes, we
work in C
symmetry. Taking the optimized geometry
for the ground state, the input is:
&STRUCTURE &END
Method casscf
End of input
&SEWARD &END
Title
1,3,-cyclopentadiene. Optimized CASSCF 431/2 geometry.
Symmetry
X
Basis set
C.ANO-L...4s3p1d.
C1 .0000000000 .0000000000 -2.3483809840
C2 .0000000000 2.2244307403 -.5647146594
C3 .0000000000 1.3954349667 1.8428732576
C4 .0000000000 -1.3954349667 1.8428732576
C5 .0000000000 -2.2244307403 -.5647146594
End of basis
Basis set
H.ANO-L...2s.
H1 1.6589565171 .0000000000 -3.5771781848
H2 .0000000000 4.1624710575 -1.1752374238
H3 .0000000000 2.5521178749 3.5144851606
H4 .0000000000 -2.5521178749 3.5144851606
H5 .0000000000 -4.1624710575 -1.1752374238
End of basis
End of Input
&SCF &END
TITLE
cyclopentadiene molecule
OCCUPIED
15 3
ITERATIONS
40
END OF INPUT
&RASSCF &END
TITLE
cyclopentadiene molecule 1A1
SYMMETRY
1
SPIN
1
NACTEL
6 0 0
INACTIVE
15 0
RAS2
0 5
THRS
1.0E-05,1.0E-03,1.0E-03
ITER
50,25
LUMORB
END OF INPUT
&ALASKA &END
NoInvariance
End of Input
&SLAPAF &END
Internal coordinates
b1 = Bond C1 C2
b2 = Bond C2 C3
b3 = Bond C3 C4
b4 = Bond C4 C5
b5 = Bond C1 H1
b6 = Bond C2 H2
b7 = Bond C3 H3
b8 = Bond C4 H4
b9 = Bond C5 H5
a1 = Angle C1 C2 C3
a2 = Angle C2 C3 C4
a3 = Angle C3 C4 C5
a4 = Angle C1 C2 H2
a5 = Angle C2 C3 H3
a6 = Angle C3 C4 H4
a7 = Angle C4 C5 H5
a8 = Angle C2 C1 H1
d1 = Dihedral H1 C1 C2 C3
Vary
b1
b2
b3
b4
b5
b6
b7
b8
b9
a1
a2
a3
a4
a5
a6
a7
a8
d1
End of Internal
Iterations
0
Thrs
0.5D-06 1.0D-03
Numerical Hessian
End of Input
We have prepared now the inputs in C symmetry and consequently
the definition of the coordinates must include more variables.
Remember that you can skip this definition and let the program
work in non-redundant coordinates. When a numerical Hessian is performed,
the keyword ITERations in SLAPAF input must be set to 0
if one does not want to proceed with a transition state search after the
calculation of the Hessian.
The Hessian obtained and the subsequent harmonic analysis are valid in a stationary point situation like here. If the Hessian is going to be used in a transition state search the keyword TS must be also included in the SLAPAF input as it will be shown in section 4.5.
If one is using AUTOMOLCAS, one must remember that the optimization
and the C Hessian cannot be computed in the same calculation.
AUTOMOLCAS will only accept successive optimization and Hessian if
both of them retain the symmetry and if the Hessian is performed
at the optimized geometry.
Cyclopentadiene in C symmetry has 18 degrees of freedom. That
means that the numerical Hessian will take 37 macro-iterations.
Finally, we obtain the force constant matrix, the harmonic
frequencies and the corresponding normal modes:
Harmonic frequencies in cm-1
1 2 3 4 5 6
850.85 856.02 970.47 1049.84 1068.32 1175.00
.08321 .20832 -.64537 -.25427 .67791 .07363
-.07705 -.01386 -.00201 -.09352 .19568 -.19069
.18007 -.01289 .32901 -.82000 .01953 -.00030
-.05590 .02564 -.01504 -.09029 -.16061 .16398
.00231 .00046 -.00347 -.00667 .01538 .00328
.00401 .00518 -.02141 -.02836 .02320 -.00542
.00763 -.00288 -.00196 -.03289 .01446 -.00277
.00807 .00006 -.00379 -.03429 -.00992 .00213
.00351 -.00634 -.02143 -.02816 -.01626 .00442
-.52971 .34451 .05712 -.11704 -.31202 .05320
.13525 -.55541 -.18763 .10262 .08994 .00053
.34009 .49441 -.22521 .08655 -.07547 .00178
.05922 -.36155 -.38406 .07864 .34060 .55013
.23422 .32449 -.10506 .27490 .22763 -.40764
-.44675 -.13619 .34222 -.35653 .28486 -.35546
.44714 -.11395 .32333 .03708 .01142 .56366
-.18748 -.10866 -.04792 -.01013 .29016 .10587
-.20940 -.06750 -.04884 .00902 .17374 .06583
7 8 9 10 11 12
1189.17 1398.34 1473.18 1494.92 1641.89 1701.33
...
13 14 15 16 17 18
1713.68 3304.88 3306.80 3317.87 3330.11 3339.58
...
There are no imaginary frequencies and therefore the geometry corresponds
to an stationary point within the C symmetry.
Sometimes may be interesting to follow the path of the optimization by looking at each one of the output files generated by STRUCTURE. The output files of each complete iteration are stored in the $WorkDir directory under the names $RANDOM.save.$iter, for instance: 2222.save.1, 2222.save.2, etc. You should not remove the $WorkDir directory if you want to keep them.
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. 4.3
in next section). The ground state has
symmetry: 1 
A. The two lowest valence excited states
are 2A and 1B
. If we optimize the geometries within
the 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 C, or even a restricted one in C, 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 calculation; it is not possible to
make state-average CASSCF calculations.
We are going to optimize the three states of thiophene in symmetry. The inputs are:
&STRUCTURE &END
Method casscf
End of input
&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
&SCF &END
TITLE
Thiophene molecule
OCCUPIED
11 1 7 3
ITERATIONS
40
END OF INPUT
&RASSCF &END
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
LUMORB
END OF INPUT
&ALASKA &END
NoInvariance
End of Input
&SLAPAF &END
Internal coordinates
b1 = Bond S1 C1
b2 = Bond C1 C2
b3 = Bond C2 C2(xy)
b4 = Bond C1 H1
b5 = Bond C2 H2
a1 = Angle S1 C1 C2
a2 = Angle S1 C1 H1
a3 = Angle C1 C2 H2
Vary
b1
b2
b3
b4
b5
a1
a2
a3
End of Internal
Iterations
20
Thrs
0.5D-06 1.0D-03
End of Input
for the ground state. For the two excited states we will replace the RASSCF inputs with
&RASSCF &END
TITLE
Thiophene molecule 2 1A1
SYMMETRY
1
SPIN
1
NACTEL
6 0 0
INACTIVE
11 0 7 1
RAS2
0 2 0 3
CIROOT
1 2
2
LEVSHFT
1.0
ITER
50,25
LUMORB
END OF INPUT
for the 2A state and
&RASSCF &END
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
LUMORB
END OF INPUT
for the 1B state.
The optimization of these three states do not present any
type of problem. But if we try to make the optimization of
the excited states or compute the numerical Hessian in
C symmetry we have convergence problems, especially
for the 1B state, which becomes now the 3A
state,
that is, the third root in the symmetry.
To help the program to converge we include one or more initial RASSCF
inputs in the input file. They will be executed sequentially
by STRUCTURE taking the final orbitals of one as starting orbitals
for the following. The final wave function in each macro-iteration
will be that of the last RASSCF. The following is an example for the calculation
of the numerical Hessian (C symmetry) of the 3A state of thiophene
at its optimized geometry.
&STRUCTURE &END
Method casscf
Do always
End of input
&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
&SCF &END
TITLE
Thiophene molecule
OCCUPIED
18 4
ITERATIONS
40
END OF INPUT
&RASSCF &END
TITLE
Thiophene molecule 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
LUMORB
END OF INPUT
&RASSCF &END
TITLE
Thiophene molecule 3 1A'
SYMMETRY
1
SPIN
1
NACTEL
6 0 0
INACTIVE
18 1
RAS2
0 5
CIROOT
1 3
3
LEVSHFT
1.0
ITER
50,25
JOBIPH
END OF INPUT
&ALASKA &END
NoInvariance
End of Input
&SLAPAF &END
Internal coordinates
b1 = Bond S1 C1
b2 = Bond C1 C2
b3 = Bond C2 C2b
b4 = Bond C2b C1b
b5 = Bond C1 H1
b6 = Bond C2 H2
b7 = Bond C1b H1b
b8 = Bond C2b H2b
a1 = Angle S1 C1 C2
a2 = Angle C1 C2 C2b
a3 = Angle C2 C2b C1b
a4 = Angle S1 C1 H1
a5 = Angle C1 C2 H2
a6 = Angle C2 C2b H2b
a7 = Angle C2b C1b H1b
Vary
b1
b2
b3
b4
b5
b6
b7
b8
a1
a2
a3
a4
a5
a6
a7
End of Internal
Iterations
0
Numerical Hessian
End of Input
Here we have two RASSCF inputs. The first one computes three state-average roots and the second one takes the orbitals generated by the first RASSCF as initial set to compute the third root. This final wave function is then used by the ALASKA and SLAPAF programs. Keyword Do always within the STRUCTURE namelist input ensures that the two RASSCF calculations will be executed in all the STRUCTURE macro-iterations. If the keyword is absent, only the last RASSCF input will be taken after the first macro-iteration.
After 31 iterations and the calculation of the force constant matrix,
the analysis of the harmonic frequencies report one imaginary value:
1144.93i. That means that although we are in a stationary point in the
symmetry, we are in a saddle point in the C hypersurface and
we must continue with the optimization in the lower symmetry.
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. It would be advantageous if excited state optimizations could be performed with state-averaged RASSCF wave functions. Codes to do this are under development.
Presently, MOLCAS is prepared to work in the point groups
C, C
, C, C, D, C
, C
, and D.
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
D
point group we have to generate the integrals and
wave function in D 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, NH
. There is
a planar transition state along the isomerization barrier between
two pyramidal structures. We want to optimize the planar structure
restricted to the D
point group. However, the electronic wave function will
be computed in C symmetry ( 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:
&STRUCTURE &END Method casscf End of input &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 &SCF &END Title NH3, planar Occupied 4 1 Iterations 40 End of Input &RASSCF &END Title NH3, planar Symmetry 1 Spin 1 Nactel 8 0 0 INACTIVE ORBITALS 1 0 RAS2 ORBITALS 6 2 LUMORB ITER 50,20 End of Input &ALASKA &END End of input &SLAPAF &END Internal coordinates bond1 = Bond N H1 bond2 = Bond N H2 bond3 = Bond N H3 ang1 = Angle H1 N H2 ang2 = Angle H1 N H3 Vary r1 = 1.0 bond1 + 1.0 bond2 + 1.0 bond3 Fix r2 = 1.0 bond1 - 1.0 bond2 r3 = 1.0 bond1 - 1.0 bond3 ang1 ang2 End of internal Iterations 20 End of input
All four atoms are in the same plane.
Working in C, 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 
equal to the sum of all three
bonds; then, we define coordinates 
and 
and keep them fixed.
will ensure that 
is equal to 
and will assure that
is equal to . and will have a zero value.
In this way all three bonds will have the same length.
As we want the system constrained into the D point group,
the three angles must be equal with a value of 120 degrees. This is
their initial value, therefore we simply keep coordinates 
and 
fixed. The result is a D structure:
******************************************** * Values of primitive internal coordinates * ******************************************** ~ BOND1 : Bond Length= 1.8957/ bohr BOND2 : Bond Length= 1.8957/ bohr BOND3 : Bond Length= 1.8957/ bohr ANG1 : Angle= 120.0000/degree, 2.0944/rad ANG2 : 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.
Another common situation in geometry optimizations is to have one or
several coordinates fixed and vary the remaining coordinates.
As an example let's 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
degrees with
respect to the other [18]. We want to optimize the
dihedral angle relating both benzene units but keep all the other
coordinates fixed. The input file should be:
&STRUCTURE &END
Method RASSCF
End of input
&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
&SCF &END
TITLE
Biphenyl twisted D2
OCCUPIED
12 9 9 11
ITERATIONS
50
END OF INPUT
&RASSCF &END
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
LUMORB
END OF INPUT
&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
Fix
b1
b2
b3
b4
h1
h2
h3
a1
a2
a3
a4
d1
d2
d3
End of Internal
Iterations
30
End of Input
To be able to optimize the molecule in that way a D 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 ( 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.
Unfortunatelly if some coordinate has to be kept frozen the
only available implementation is to use the user-defined
internal coordinates. In future releases it will be possible
to froze coordinates in the non-redundant internal approach.
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 (
180 for Dihedral angles and
90 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. Taking the
previous example on biphenyl, when we try to use the geometry of
the input and vary also bond distances, the program respond with
the message:
Diverging in INTOCA!!!
4 B4 -.8227E+00
rMax= -0.822686317639044873
rOld= -0.171729836976677319E-01
This is an indication of a problem in the primitive internal coordinate B4, in this case a C-C bond. Options QUASi trust radius or MAXStep may help to solve it, although it is probably better to start from a better geometry. 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.
