8
$\begingroup$

I am computing the optimized molecular structure of an organic molecule in the Gaussian 09 computational chemistry package.
As I have it now, I do not think that the optimized structure that I am getting from Gaussian is correct, in part due to the symmetry of some of the bonds. To try and correct this, I want to fix two of the bond angles in the molecule in order to help/force it to produce a more believable structure. Even though I want to fix these bond angles, I want to leave everything else variable and open to optimization, including those bond lengths.

How can I specify this in the input file?

Fixing the bond angle to a set value seems like an extreme solution to this problem. Perhaps I can just constrain the angle within a certain range, is this possible?

After a lot of work and playing with settings, I have figured out how to fix bond angles, but now my difficulty is fixing them to 180 degrees (my desired value), for linear bonds. It seems that this is very hard to do.

Fixing the bond angle to 179 degrees (or 179.9 degrees), results in an error saying "Error in internal coordinate system." A lot of people online say this is because of linear/nearly-linear bonds in the molecule being optimized. Some people also say this might be avoided if I try to optimize in Cartesian coordinates, instead of internal coordinates, using the cartesian parameter after opt. However, I do not think I can use ModRedundant and cartesian parameters at the same time.

How can I fix a bond angle to be linear, but still optimize the rest of the structure in Cartesian coordinates? Is it even possible to fix angles to 180 degrees?

Improve this question
$\endgroup$
3
  • $\begingroup$ Is this an organic or inorganic structure and what method do you use? But fixing an angle to sth that you believe might not be what is "right". $\endgroup$
    – pH13 - Yet another Philipp
    Commented Oct 15, 2015 at 7:05
  • $\begingroup$ And what method do you use for that organic molecule, which leads to conpletely unbelievable results? $\endgroup$
    – pH13 - Yet another Philipp
    Commented Oct 17, 2015 at 2:29
  • $\begingroup$ I don't think it is the method/basis set's fault. Try many different choices all give me the same problem. $\endgroup$
    – user38694
    Commented Oct 19, 2015 at 0:16

4 Answers 4

Reset to default
9
$\begingroup$

Fixing the bond angle to a set value seems like an extreme solution to this problem.

It is. To be honest, fixing parameters in a calculation is almost always no solution at all. The resulting geometry has little to no meaning and all derived properties are for a purely hypothetical case. If the located geometry is no stationary point, what does that even mean? A wrong solution is worse than no solution at all. You have to be very careful with what you fix and when you fix, and you should always confirm what you have found with a full optimisation and a frequency calculation. (On various levels of theory.)

I think Greg's advice is also very important. Symmetry can be a very powerful tool if utilised correctly, but you should make sure that the implied symmetry is correct. You also need to check whether your methodology matches the system you are computing.

All that being said, there are certainly a few very productive uses for partial optimisations:

  • finding higher or unusual local minima
  • scanning the potential energy surface
  • modifying converged geometries
  • fixing transitional modes
  • etc.

There are a few possibilities how to make partial optimisations in gaussian. All have some caveats and all may lead to undesired results. There are certainly restrictions to all of them. I'll give a short overview the most common ones with a few practical examples.

You have already noticed, that fixing a bond angle to 180 degrees is a challenge. This is due to the fact, that it is a requirement to have linearly independent variables, which is simply not possible with threes point on a line. However, there are certain ways to work around the problem. Spoiler alert: It is only possible with z-matrices or very rigid symmetry (where you fix the symmetry, not the bond angle).

Partial optimisations with z-matrices

Let's start with something very simple. A full optimisation of methanol without any symmetry restrictions.

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
opt
int(ultrafinegrid)
sym(none)

title card required

0 1
 c
 h    1 hc2
 h    1 hc3         2 hch3
 h    1 hc4         2 hch4          3 dih4
 o    1 oc5         2 och5          3 dih5
 h    5 ho6         1 hoc6          2 dih6

hc2           1.1
hc3           1.1
hch3        109.5
hc4           1.1
hch4        109.5
dih4        120.0
oc5           1.4
och5        109.5
dih5       -120.0
ho6           0.9
hoc6        109.5
dih6        180.0

This will result in the following optimised geometry:

 c
 o    1 oc2
 h    1 hc3         2 hco3
 h    1 hc4         2 hco4          3 dih4
 h    1 hc5         2 hco5          4 dih5
 h    2 ho6         1 hoc6          3 dih6

oc2           1.414851
hc3           1.109697
hco3        107.476
hc4           1.118853
hco4        113.586
dih4       -118.329
hc5           1.118853
hco5        113.586
dih5       -123.342
ho6           0.974680
hoc6        107.408
dih6        179.999

Fixing an angle in a z-matrix is straight forward. You need to specify that you are doing a partial optimisation (popt) and move the fixed parameter to the constants section. As an example, I am fixing the COH angle, and also invoke symmetry.

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
popt
int(ultrafinegrid)

title card required

0 1
 c
 h    1 hc2
 h    1 hc3         2 hch3
 h    1 hc3         2 hch3          3 dih4
 o    1 oc5         2 och5          3 dih5
 h    5 ho6         1 coh           2 dih6
  Variables:
hc2           1.100000
hc3           1.100000
hch3        109.5
dih4        120.000
oc5           1.400000
och5        109.5
dih5       -120.000
ho6           0.900000
dih6        180.000
  Constants:
coh          90.0

This will give you the following list of optimised parameters:

       hc2         1.1104 
       hc3         1.1173 
      hch3       107.0973 
      dih4       116.0195 
       oc5         1.4471 
      och5       108.92   
      dih5      -121.9694 
       ho6         0.9796 
      dih6       180.0092 
       coh        90.0

Now we will try the same, but we will make the bond angle 180 degrees. Here we need to introduce dummy atoms:

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
popt
int(ultrafinegrid)

title card required

0 1
 c
 h    1 hc2
 h    1 hc3         2 hch3
 h    1 hc3         2 hch3          3 dih4
 o    1 oc5         2 och5          3 dih5
xx    5 xxo6        1 xxoc6         2 dih6
 h    5 ho7         6 xxoc6         1 dih7
  Variables
hc2           1.100000
hc3           1.100000
hch3        109.500
dih4        120.000
oc5           1.400000
och5        109.500
dih5       -120.000
dih6        180.000
ho7           0.950000
  Constants
xxo6          1.000000
xxoc6        90.000
dih7        180.000

Which will give you once again the following set of parameters:

       hc2         1.1209
       hc3         1.1211
      hch3       106.3981
      dih4       113.09  
       oc5         1.3644
      och5       112.3772
      dih5      -123.4351
      dih6       180.0   
       ho7         0.9465
      xxo6         1.0   
      xxoc6       90.0   
      dih7       180.0

And the following Cartesian coordinates:

    6
scf done:  -115.581115
 C     0.000000     0.000000     0.000000
 H     0.000000     0.000000     1.120932
 H     1.075500     0.000000    -0.316498
 H    -0.421786    -0.989342    -0.316498
 O    -0.695163     1.052866    -0.519431
 H    -1.177401     1.783245    -0.879763

Quick and Dirty fixing with Cartesian coordinates

Appending -1 to an element in a row in the Cartesian coordinates block will fix it. This is a very quick and dirty approach, which may lead to a brute force computation. You are basically restricting the algorithm to rotate the molecule, which leads to less effective calculations. Fixing Cartesian coordinates also fixes the bond length, but it is the only way to force a 180 degree angle with this choice of input. Here is an example input:

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
popt
int(ultrafinegrid)
symm(loose)

title card required

0 1
 C  -1    0.000000     0.000000     0.000000
 H        0.000000     0.000000     1.300000
 H        1.221600     0.000000    -0.444626
 H       -0.479083    -1.123738    -0.444626
 O  -1   -0.695162     1.052869    -0.519426
 H  -1   -1.177400     1.783250    -0.879755

This will result in the following optimised geometry:

    6
scf done:  -115.581115
 C    -0.663620    -0.000005     0.000000
 H    -1.090591     0.003267    -1.036538
 H    -1.090556     0.896080     0.521076
 H    -1.090648    -0.899281     0.515456
 O     0.700780    -0.000004     0.000001
 H     1.647272    -0.000004     0.000001

Note that this calculation takes about five times longer.

Fixing parameters with modredundant

This option is described quite straight forward in the Gaussian online manual. But let's make an example where we try to restrict the COH angle to 90 degrees.

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
opt(modredundant)
int(ultrafinegrid)

title card required

0 1
 C     0.000000     0.000000     0.000000
 O    -0.659966     1.143095    -0.466667
 H    -0.000000     0.000000     1.089000
 H     1.026720     0.000002    -0.362997
 H    -0.513360    -0.889165    -0.363000
 H    -0.502133     0.869719    -1.359507

A 5 6 1 F

With this option it is impossible to fix the angle to be 180 degrees, since this is not an allowed value. You also cannot use dummy atoms, as these do not exist in the redundant coordinate system. It might be a good idea, to kill this variable completely, to ensure that the program will not fail trying to construct it, or when an optimisation is getting close to that value.

Alternatively the angle can be restricted via a linear bend. This is the best option, if you know that the angle is 180.0 degrees.

%chk=bp86svp.chk
#p BP86/def2SVP/W06
DenFit
opt(Modredundant)
int(ultrafinegrid)
sym(none)

title card required

0 1
C        0.000000000      0.000000000      0.736983516
H        1.037151577      0.000000000      1.103476471
H       -0.518575789     -0.898199613      1.103476471
H       -0.518575789      0.898199613      1.103476471
O        0.000000000      0.000000000     -0.663016475
H        0.000000000      0.000000000     -1.563016262

L 1 5 6 2 F

The last line freezes the linear angle between atoms 1, 5, and 6. Atom 2 needs to be specified to determine the plane in which the linear bend lies, or is orthogonal to.

Fixing linear angles with general internal coordinates (GIC)

A new way of handling the coordinate system was introduced with Gaussian 16: generalised internal coordinates. It is in many ways superior to earlier versions, however, also a little bit harder to learn. Nevertheless, it is absolutely worthwhile trying.

Here is a sample input for oxoketene $\ce{O=C=C=O}$. In this case I specified the complete set of coordinates necessary.

#p PM6
opt geom(AddGIC)
scf(xqc,MaxCycle=250)
int(ultrafinegrid)
symmetry(none)

Oxoketene (constrained)

0 1
O        0.000000000      0.000000000     -2.210000000
C        0.000000000      0.000000000     -0.770000000
C        0.000000000      0.000000000      0.770000000
O        0.000000000      0.000000000      2.210000000

BondC1O1(Active)=R(2,1)
BondC2O2(Active)=R(3,4)
BondC1C2(Active)=R(2,3)
LinBend1O1C1C2(Freeze)=L(1,2,3,0,-1)
LinBend2O1C1C2(Freeze)=L(1,2,3,0,-2)
LinBend1C1C2O2(Freeze)=L(2,3,4,0,-1)
LinBend2C1C2O2(Freeze)=L(2,3,4,0,-2)

The format is 

L(     i,      j,      k,      l,                                 M)
L(Atom 1, Atom 2, Atom 3, <Atom 4 | Orthogonal to Plane>, component)

The linear bend is therefore defined between i, j, k as Atom 1, Atom 2, Atom 3 and must not be dummy atoms (and with GIC these are hardly necessary anymore). The fourth value l, can be another atom to determine the orthogonality/ plane to, or it can be defined as a plane where l = 0, -1, -2, -3 standing for automatic, $yz$, $xz$, $xy$, respectively. The last value M gives the component, of which two are necessary, and must be M = -1, -2.

The methanol example is as follows:

#p PM6
DenFit
opt(AddGIC)
int(ultrafinegrid)
sym(none)

title card required

0 1
C        0.000000000      0.000000000      0.736983516
H        1.037151577      0.000000000      1.103476471
H       -0.518575789     -0.898199613      1.103476471
H       -0.518575789      0.898199613      1.103476471
O        0.000000000      0.000000000     -0.663016475
H        0.000000000      0.000000000     -1.563016262

L(1, 5, 6, 2, -1) Freeze
L(1, 5, 6, 2, -2) Freeze
Improve this answer
$\endgroup$
4
$\begingroup$

As mentioned, Z-matrix input is one solution. You also can constrain angle (and other vairable) of Cartesian structures using opt(ModRedundant) Gaussian manual

If your problem is related to symmetry, you should be careful about two other things:

  • Do you enforce a (probably incorrect) symmetry? Gaussian may picks up a wrong symmetry due to a bad input geometry. You can check that from your output file, and see what symmetry are you constrained to. You can switch off symmetry (both for optimization and SCF part) or can modify your input that Gaussian can find the right symmetry.

  • Do you have the right method to represent your wavefunction? It is often found mistake that you you a single determinant method to a multideterminant problem or you don't realize the presence of static or dynamic Jahn-Teller effect. If your structure is very different from the expected, you may have this problem.

About related question: Not that I am aware off, but generally it makes not much sense anyway. If you want to fix e.g. the HOH angle in water between 80 and 90 deg, it will be 90 anyways - same as you fix it at 90 deg. I.e. it will stick to one extreme, that you already can figure out from chemical intuition or previous, unconstrained runs.

Improve this answer
$\endgroup$
3
$\begingroup$

Thanks for all answers it really helped me. By the way, all links which has given above were died so here is the new gaussian manual link for z-matrix: http://gaussian.com/zmat/

I generated my z-matrix file by using Avagadro (from extensions->> Gaussian menu) software here is the download link: https://sourceforge.net/projects/avogadro/files/latest/download I hope it will work for somebody.

Improve this answer
$\endgroup$
1
$\begingroup$

I believe one can accomplish this by giving the geometry input in z-matrix form, and using the 'constants' and 'variables' keywords. This page from the Gaussian 09 manual shows in detail how one could use these keywords.

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.