7
$\begingroup$

Sometimes it is necessary to perform a partial optimisation of a molecule, i.e. keeping certain variables constant. One example might be to pre-optimise a transition state, taking advantage of not having to compute the force constants in the initial set-up.

For example, let's consider the transition state of the Diels-Alder Reaction of ethene and butadiene, where we want to freeze the distance between the two molecules.

How can I achieve this in the software suite GAMESS?

$\endgroup$
10
$\begingroup$

Unfortunately GAMESS is not as user friendly as other quantum chemistry packages, so when it comes to running non-standard calculations a little bit more effort is necessary.

First of all I recommend optimising the reactants by themselves, so that you have a somewhat valid idea what the bond lengths on your level of theory are going to be. After that you should assemble your best guess for the structure. In this example I have chosen a separation between the molecules of about 220 pm.

Diels-Alder TS guess

The coordinates of this structure are:

! DIENE
C       -0.486294808     -0.431915730     -1.395523790
C       -1.166676574      0.550508040     -0.705577944
C       -1.166676574      0.550508040      0.705577944
C       -0.486294808     -0.431915730      1.395523790
H       -0.352988862     -0.347708199      2.465648783
H       -0.084717264     -1.289523730      0.864198908
H       -1.712490163      1.338624758      1.223514576
H       -1.712490163      1.338624758     -1.223514576
H       -0.352988862     -0.347708199     -2.465648783
H       -0.084717264     -1.289523730     -0.864198908
! DIENOPHILE
C        1.559305232     -0.085874351     -0.690651129
C        1.559305232     -0.085874351      0.690651129
H        1.514911897     -1.022640981      1.232100569
H        1.603704743      0.851022589      1.232175888
H        1.603704743      0.851022589     -1.232175888
H        1.514911897     -1.022640981     -1.232100569

There are approaches, where you would simply fix the Cartesian coordinates. Here we only want to fix two bond lengths, 1(C)-11(C) and 4(C)-12(C). If we were to fix the coordinates of these atoms, we would also fix 1(C)-4(C) and 11(C)-12(C). Therefore a better approach is to use internal coordinates, in GAMESS the keyword group $ZMAT is used.
Since there are a lot, it is tricky to set them up. We will therefore use the implemented automated procedure. Starting with a dry run, i.e. EXETYP=CHECK, so that we can check what the program actually does. We need to define $3\times N-6=42$ internal coordinates, thus NZVAR=42. Also note, that we do not use symmetry in this case. Because of the nature of the structure we need to add some bonds already manually, NONVDW(1)=1,11, 4,12. This is the resulting input:

 $CONTRL
   SCFTYP=RHF
   RUNTYP=ENERGY
   EXETYP=CHECK
   COORD=UNIQUE
   ICHARG=0
   MULT=1
   NZVAR=42
   $END
 $BASIS
   GBASIS=STO
   NGAUSS=3
   $END
 $ZMAT
   AUTO=.TRUE.
   DLC=.TRUE.
   NONVDW(1)=1,11, 4,12
   $END
 $GUESS
   GUESS=HUCKEL
   $END
 $DATA
   C6H10
   C1
   CARBON      6.0     -0.486294808        -0.431915730        -1.395523790
   CARBON      6.0     -1.166676574         0.550508040        -0.705577944
   CARBON      6.0     -1.166676574         0.550508040         0.705577944
   CARBON      6.0     -0.486294808        -0.431915730         1.395523790
   HYDROGEN    1.0     -0.352988862        -0.347708199         2.465648783
   HYDROGEN    1.0     -0.084717264        -1.289523730         0.864198908
   HYDROGEN    1.0     -1.712490163         1.338624758         1.223514576
   HYDROGEN    1.0     -1.712490163         1.338624758        -1.223514576
   HYDROGEN    1.0     -0.352988862        -0.347708199        -2.465648783
   HYDROGEN    1.0     -0.084717264        -1.289523730        -0.864198908
   CARBON      6.0      1.559305232        -0.085874351        -0.690651129
   CARBON      6.0      1.559305232        -0.085874351         0.690651129
   HYDROGEN    1.0      1.514911897        -1.022640981         1.232100569
   HYDROGEN    1.0      1.603704743         0.851022589         1.232175888
   HYDROGEN    1.0      1.603704743         0.851022589        -1.232175888
   HYDROGEN    1.0      1.514911897        -1.022640981        -1.232100569
   $END

The resulting auxiliary outputfile .dat has now the group that we will use from here on in. And the .log file will give us the exact value of the internal coordinate which we want to freeze.

                     --------------------
                     INTERNAL COORDINATES
                     --------------------

                 - - ATOMS - -         COORDINATE      COORDINATE
 NO.   TYPE    I  J  K  L  M  N        (BOHR,RAD)       (ANG,DEG)
 ----------------------------------------------------------------
   1 STRETCH   1 11                     4.1406435       2.1911343
...
  10 STRETCH   4 12                     4.1406435       2.1911343
...

We can now modify our $ZMAT group accordingly. Unfortunately the input gets very long, if you would like to retain some readability. Note that we create redundant internal coordinates, so we need to increase the number of NZVAR, too, to read all of them in.
Technically we do not have to set the generation of coordinates to false, but it helps remembering what you have done if you come back in a month or so, therefore AUTO=.FALSE.. Now we add the coordinates we would like to freeze with IFZMAT(1)=<list> and the corresponding values FVALUE(1)=<list>.

 $CONTRL 
   SCFTYP=RHF 
   RUNTYP=OPTIMIZE
   COORD=UNIQUE
   ICHARG=0
   MULT=1
   NZVAR=89  
   $END
 $BASIS  
   GBASIS=STO 
   NGAUSS=3
   $END
 $ZMAT
   AUTO=.FALSE.
   DLC=.TRUE.
   IZMAT(1)=
   1,   1,  11,
   1,   1,  10,
   1,   1,   9,
   1,   1,   2,
   1,   2,   8,
   1,   2,   3,
   1,   3,   7,
   1,   3,   4,
   1,   4,   6,
   1,   4,  12,
   1,   4,   5,
   1,  11,  16,
   1,  11,  15,
   1,  11,  12,
   1,  12,  14,
   1,  12,  13,
   2,   1,  11,  15,
   2,   1,  11,  16,
   2,   1,   2,   3,
   2,   1,   2,   8,
   2,   1,  11,  12,
   2,   2,   3,   7,
   2,   2,   1,  10,
   2,   2,   1,  11,
   2,   2,   1,   9,
   2,   2,   3,   4,
   2,   3,   2,   8,
   2,   3,   4,   5,
   2,   3,   4,  12,
   2,   3,   4,   6,
   2,   4,   3,   7,
   2,   4,  12,  13,
   2,   4,  12,  14,
   2,   4,  12,  11,
   2,   5,   4,  12,
   2,   5,   4,   6,
   2,   6,   4,  12,
   2,   9,   1,  10,
   2,   9,   1,  11,
   2,  10,   1,  11,
   2,  11,  12,  13,
   2,  11,  12,  14,
   2,  12,  11,  15,
   2,  12,  11,  16,
   2,  13,  12,  14,
   2,  15,  11,  16,
   3,   1,   2,   3,   4,
   3,   1,   2,   3,   7,
   3,   1,  11,  12,  14,
   3,   1,  11,  12,   4,
   3,   1,  11,  12,  13,
   3,   2,   1,  11,  15,
   3,   2,   1,  11,  16,
   3,   2,   3,   4,   5,
   3,   2,   3,   4,  12,
   3,   2,   3,   4,   6,
   3,   2,   1,  11,  12,
   3,   3,   2,   1,  11,
   3,   3,   4,  12,  13,
   3,   3,   4,  12,  14,
   3,   3,   4,  12,  11,
   3,   3,   2,   1,   9,
   3,   3,   2,   1,  10,
   3,   4,  12,  11,  15,
   3,   4,  12,  11,  16,
   3,   4,   3,   2,   8,
   3,   5,   4,  12,  11,
   3,   5,   4,  12,  14,
   3,   5,   4,  12,  13,
   3,   5,   4,   3,   7,
   3,   6,   4,  12,  11,
   3,   6,   4,  12,  14,
   3,   6,   4,  12,  13,
   3,   6,   4,   3,   7,
   3,   7,   3,   2,   8,
   3,   7,   3,   4,  12,
   3,   8,   2,   1,   9,
   3,   8,   2,   1,  10,
   3,   8,   2,   1,  11,
   3,   9,   1,  11,  16,
   3,   9,   1,  11,  15,
   3,   9,   1,  11,  12,
   3,  10,   1,  11,  12,
   3,  10,   1,  11,  15,
   3,  10,   1,  11,  16,
   3,  13,  12,  11,  15,
   3,  13,  12,  11,  16,
   3,  14,  12,  11,  15,
   3,  14,  12,  11,  16,
   IFZMAT(1)=
   1,   1,  11,
   1,   4,  12,
   FVALUE(1)=
   2.1911343,
   2.1911343,
   $END
 $GUESS  
   GUESS=HUCKEL
   $END
 $DATA
   C6H10
   C1
   CARBON      6.0     -0.486294808        -0.431915730        -1.395523790
   CARBON      6.0     -1.166676574         0.550508040        -0.705577944
   CARBON      6.0     -1.166676574         0.550508040         0.705577944
   CARBON      6.0     -0.486294808        -0.431915730         1.395523790
   HYDROGEN    1.0     -0.352988862        -0.347708199         2.465648783
   HYDROGEN    1.0     -0.084717264        -1.289523730         0.864198908
   HYDROGEN    1.0     -1.712490163         1.338624758         1.223514576
   HYDROGEN    1.0     -1.712490163         1.338624758        -1.223514576
   HYDROGEN    1.0     -0.352988862        -0.347708199        -2.465648783
   HYDROGEN    1.0     -0.084717264        -1.289523730        -0.864198908
   CARBON      6.0      1.559305232        -0.085874351        -0.690651129
   CARBON      6.0      1.559305232        -0.085874351         0.690651129
   HYDROGEN    1.0      1.514911897        -1.022640981         1.232100569
   HYDROGEN    1.0      1.603704743         0.851022589         1.232175888
   HYDROGEN    1.0      1.603704743         0.851022589        -1.232175888
   HYDROGEN    1.0      1.514911897        -1.022640981        -1.232100569
   $END

The optimisation should look something like the following animation; and took 19 steps in my example.

animation of optimisation

Now you could use the results of this calculation to run a full fledged transition state search. Remember, that you have all necessary information to restart a run in the auxiliary output, i.e. the .dat file. You should always keep this one until you are completely done with your calculations.

If you want to calculate along and wonder what the optimised geometry is, here are the results of the final step:

----- RESULTS FROM SUCCESSFUL RHF      GEOMETRY SEARCH -----
----- COORDS, ORBS, GRADIENT, AND APPROX. HESSIAN -----
 COORDINATES OF SYMMETRY UNIQUE ATOMS (ANGS)
   ATOM   CHARGE       X              Y              Z
 ------------------------------------------------------------
 CARBON      6.0  -0.4847330011  -0.4122786646  -1.4086312155
 CARBON      6.0  -1.1902020846   0.5288839653  -0.7019679480
 CARBON      6.0  -1.1902020983   0.5288839622   0.7019679173
 CARBON      6.0  -0.4847329728  -0.4122786820   1.4086312270
 HYDROGEN    1.0  -0.4076405265  -0.3274824211   2.4860530248
 HYDROGEN    1.0  -0.4019303863  -1.4244097276   1.0396872383
 HYDROGEN    1.0  -1.5467259917   1.4139143889   1.2140730563
 HYDROGEN    1.0  -1.5467259527   1.4139143906  -1.2140730901
 HYDROGEN    1.0  -0.4076405381  -0.3274824126  -2.4860529903
 HYDROGEN    1.0  -0.4019304494  -1.4244096950  -1.0396872193
 CARBON      6.0   1.5587691569  -0.0858507824  -0.6884831136
 CARBON      6.0   1.5587691834  -0.0858507803   0.6884831288
 HYDROGEN    1.0   1.9167366148  -0.9561957837   1.2206602417
 HYDROGEN    1.0   1.6758814772   0.8473264236   1.2212666939
 HYDROGEN    1.0   1.6758814263   0.8473264124  -1.2212666800
 HYDROGEN    1.0   1.9167365785  -0.9561957779  -1.2206602275
--- OPTIMIZED RHF      MO-S --- GENERATED AT Thu Dec 14 23:42:52 2017
E=     -230.0376739583, E(NUC)=  227.2020852580
 $VEC
...

Disclaimer: There is another way of freezing coordinates. I might add that at another point as another answer. It is a weeny bit more complicated, however offers also some more flexibility.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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