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?
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.
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.
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.