7
$\begingroup$

I would like to perform an optimization in Gaussian 09 of my molecule with two dihedral angles being frozen. I would like that the values of these dihedral angles are not fixed, but can relax within a range of ±5°. In Gaussian 03 it seems to be possible using Opt=ModRedundant keyword and the following syntax:

[Type] N1 [N2 [N3 [N4]]] [[+=]value] [A | F] [[min] max]]

I have tried using the same syntax in Gaussian 09, but it doesn't work. Is it possible to perform such calculation in Gaussian 09 or 16 without using the scan approach?

The calculation fails with the following error:

The following ModRedundant input section has been read: 
Invalid extra data found.
$\endgroup$
1
  • $\begingroup$ Do you know it is possible in G03 or do you just think so? $\endgroup$ Commented Sep 21, 2018 at 17:07

2 Answers 2

5
$\begingroup$

TL;DR: It is impossible (=not documented if it were) with Gaussian 09 Rev. E.01 and older, it might be possible with Gaussian 16, but in all cases not via opt(modredundant).

I am unaware of any sloppy freezing possibility in pre-Gaussian 16, and for the latter one has to dig much deeper in the manual and start constructing it; but at least it works in theory.

Gaussian 03

Unfortunately, you have misinterpreted the commands you found in the manual. You can find the full description for the opt keyword in via The Internet Archive, but I'll share the important parts here for convenience.

The construct

[Type] N1 [N2 [N3 [N4]]] [[+=]value] [A | F] [[min] max]]

indeed exists for a opt(modredundant) input section, but it behaves not as a sloppy freeze, but as a range exclusion criterion, which you can find a little lower in the manual:

[...]
An asterisk (*) in the place of an atom number indicates a wildcard. Min and max then define a range (or maximum value if min is not given) for coordinate specifications containing wildcards. The action specified by the action code is taken only if the value of the coordinate is in the range.
[...]

For example, the input

D * * * * F 175.0

freezes all defined dihedral angles that are smaller than 175.0°.

The input

D * 4  5 * F 160.0 175.0

freezes all defined dihedral angles around the bond between atom 4 and 5 in the interval 160.0 - 175.0°.

The input

D 3 4  5 6 F 160.0 175.0

is a syntax error because it does not contain a wildcard.

Gaussian 09

The option modredundant to the opt keyword still exists, but has a slightly modified input. Specifically the [[+=]value]] part has been removed. Again you can look up the full page on The Internet Archive.
The examples above should, however, still behave the same.

[...]
Lines in a ModRedundant input section use the following syntax:
[Type] N1 [N2 [N3 [N4]]] [A | F] [[min] max]]
[...]
An asterisk (*) in the place of an atom number indicates a wildcard. Min and max only apply to coordinate specifications containing wildcards. The action then specified by the action code is taken only if the value of the coordinate is in the minmax range (or below maximum value if min is not given).
[...]

Therefore you also still get the error "Invalid extra data found".

Gaussian 16

In this version, while opt(modredundant) is still maintained, a new feature is rolled out. As far as I can see, there are no changes to Gaussian 09 regarding the modredundant input section. You can find it live on the Gaussian homepage or a copy on The Internet Archive.

According to this, the above examples should still produce the same results.

The new feature is Generalized Internal Coordinates or GIC. See a description on the Gaussian homepage or on The Internet Archive. You can see the GIC in action in my answer on Constrain bond angle In Gaussian molecular structure optimization and on Gaussian: Relaxed scan with modredundant optimization and dummy atoms.

In theory you should be able to freeze a GIC if it meets certain criteria for example from the manual:

The following example removes an angle coordinate generated by default if ≥179.9°, substituting a linear bend:

A(1,2,3) Remove Min=179.9        !  Remove angle coordinate if too large.  
L(1,2,3,0,-1) Add IfNot A(1,2,3) !  Add linear bends only if 
L(1,2,3,0,-2) Add IfNot A(1,2,3) !+ the angle coordinate not active.  

I have prepared a small example, but I am really not sure if it worked as intended. For cis-butadien the dihedral distortions should be small, so there isn't much visible. There are no syntax errors, but I am not really sure whether the section actually did what it should. Try it yourself and investigate further. (! recognises comments)

#p PM6                ! Semiempirical method for fast testing
opt                   ! Optimise
geom(AddGIC)          ! Use GIC, can also be specified with opt(AddGIC)
scf(xqc,MaxCycle=250) ! Switch to quadratic conv. if necessary, more cycles
symmetry(none)        ! Do not use symmetry

distorted cis-Butadiene (semi-constrained)
optimisation run
[Free title card can have multiple lines]

0 1
 C     0.000000     0.000000     0.000000
 H     0.000000     0.000000     1.089000
 H     1.026720     0.000000    -0.362996
 H    -0.513358    -0.889166    -0.363001
 C    -0.754249     1.306394    -0.533333
 H    -0.738680     2.222526     0.055206
 C    -1.531373     1.266770    -1.931370
 H    -2.044734     2.155933    -2.294370
 C    -1.615181    -0.098364    -2.761662
 H    -2.540288    -0.117779    -3.335874
 H    -0.765588    -0.163525    -3.439797
 H    -1.596709    -0.942933    -2.074433

sloppyfix=D(1,5,7,9)               ! Define an alias for the dihedral
sloppyfix freeze min=5.0 max=355.0 ! freeze in intervall

! Empty line above is important.

That the setup somewhat works can be seen with the following example of dihydrogen. First the atoms are separated 400 pm (=4.00 Å), and the value should be frozen at 250 pm (=2.50 Å), the PM6 optimised value is 76 pm (=0.76 Å).

#p PM6 opt(AddGIC,MaxStep=5) scf(xqc,MaxCycle=250) symmetry(none)

Free title card

0 1
 H    -2.000000     0.000000     0.000000
 H     2.000000     0.000000     0.000000

theHHbond=R(1,2)
theHHbond freeze max=2.5

! Empty line above is important.

Note that it will be frozen after the condition triggers, so it will not be able to the value that triggers it. Here are the values and messages you'll find in the optimisation run:

                           ----------------------------
                           !    Initial Parameters    !
                           ! (Angstroms and Degrees)  !
 --------------------------                            --------------------------
 ! Name      Definition          Value          Derivative Info.                !
 --------------------------------------------------------------------------------
 ! theHHbond R(1,2)              4.0            estimate D2E/DX2                !
 --------------------------------------------------------------------------------

[...]

 Variable       Old X    -DE/DX   Delta X   Delta X   Delta X     New X
                                 (Linear)    (Quad)   (Total)
 theHHbon     4.97283  -0.03796   0.00000  -1.00000  -1.00000   3.97283
         Item               Value     Threshold  Converged?
 Maximum Force            0.037962     0.000450     NO
 RMS     Force            0.037962     0.000300     NO
 Maximum Displacement     0.500000     0.001800     NO
 RMS     Displacement     0.707107     0.001200     NO
 Predicted change in Energy=-4.653046D-02
 Lowest energy point so far.  Saving SCF results.
 Internal coordinate with a value of  3.9728 a.u. has been conditionally frozen
 because it is in Min/Max range.
 The coordinate's label:
 theHHbond
 GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

[...]

 Variable       Old X    -DE/DX   Delta X   Delta X   Delta X     New X
                                 (Linear)    (Quad)   (Total)
 theHHbon     3.97283  -0.07253   0.00000   0.00000   0.00000   3.97283
         Item               Value     Threshold  Converged?
 Maximum Force            0.000000     0.000450     YES
 RMS     Force            0.000000     0.000300     YES
 Maximum Displacement     0.000000     0.001800     YES
 RMS     Displacement     0.000000     0.001200     YES
 Predicted change in Energy=-0.000000D+00
 Optimization completed.
    -- Stationary point found.
                           ----------------------------
                           !   Optimized Parameters   !
                           ! (Angstroms and Degrees)  !
 --------------------------                            --------------------------
 ! Name      Definition          Value          Derivative Info.                !
 --------------------------------------------------------------------------------
 ! theHHbond R(1,2)              2.1023         -DE/DX =   -0.0725              !
 --------------------------------------------------------------------------------
 Lowest energy point so far.  Saving SCF results.
 Largest change from initial coordinates is atom    1       0.949 Angstoms.
 GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

Note that in this example there is only one coordinate, ergo the optimisation finishes when the freeze condition is triggered. In this case there is quite a significant deviation from the condition max=2.5 and the actual value when it is frozen 2.1023. I would expect it to behave better with larger systems.

I would assume that the variable remains frozen for the remainder of the optimisation, with no more changes possible.
Standard disclaimer: Geometry optimisations with frozen parameters are critical for production purposes. Always check if you have a well converged molecular structure (stationary point) with no imaginary modes.

$\endgroup$
1
$\begingroup$

I have a somewhat acceptable solution to this problem. As Martin said in the other answer, it's pretty much impossible to do it with sloppy restraints imposed from the beginning, because the way Freeze works with GIC's Min and Max conditions is such. If you happen to start with the coordinate in the range you want to sloppy freeze in, Gaussian will immediately freeze it at step #1. Annoying. Instead, try this: Let's say you have an angle that you want to maintain within 0-10 degrees, and you happen to start at step #0 with the angle being 6.1 degrees. Like explained above, if you say

A(1,2,3) \ Freeze\ Min=0.0 \ Max=10.0

you will get nowhere beyond 6.1 degrees. Instead, find out which way it is headed after a few steps. Let's say the optimization goes in the way of Angle = 6, 6.5, 6.8... and so on, ending in 13.5 degrees. Change your restraints to

A(1,2,3) \ Freeze \ Min=9.6 \ Max=10.1

Now, your optimization will go in the way of Angle = 6, 6.5, 6.8.....7.6,... 9.7 and stop right there. Your angle is now frozen at 9.7 degrees. Et voila!

The real problem begins when you have a shitty system that goes 6.1, 6.5, 7.6, 5.4, 5.2, 1.0, -5, -6, ... That hasn't happened to be yet, and I'm hoping I can catch rare occurrences like these ad hoc.

New contributor
Adharsh Raghavan is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$

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