The problem: Given a starting geometry and an interaction model (i.e., a potential energy surface - PES) find the system geometry that corresponds to the lowest energy.
Bear in mind the following:
- This will not always result in the global minimum.
- Most likely you will end up with a local minimum or even a saddle point.
- A frequency analysis can tell you which of the two (minimum or saddle point) you've got.
To find a global minimum you will need to do a large number of minimizations from suitably chosen random starting points. Or, use the basin-hopping algorithm.
- The geometries of the individual interacting parts of the system. These will typically be the molecules that will be held rigid subsequently.
- The starting geometry of the whole system. For example, the positions/rotations for the interacting molecules.
- A PES: the energy model that describes how the molecules interact.
In the example used here you will need to use the data for Model(1) of the pyridine dimer. Please pick up the data from that link. Make sure that the example provided works, and that the output agrees with the sample output provided.
Orient file for geometry optimization
pyr2_geom_opt.ornt> UNITS BOHR Parameters Sites 100 polarizable 100 S-functions 50000 Alphas 50000 Parameter-sets 50000 Pairs 100000 Molecules 10 End Types H1 Z 1 H2 Z 1 H3 Z 1 H4 Z 1 H5 Z 1 N Z 7 C1 Z 6 C2 Z 6 C3 Z 6 C4 Z 6 C5 Z 6 End ! These variables define the starting geometry for the pyridine dimer ! You may alter this. Variables x1 0.052635 B y1 0.245807 B z1 -5.459405 B psi1 96.758 D nx1 0.926446 ny1 -0.277043 nz1 0.254842 x2 -0.052635 B y2 -0.245807 B z2 5.459405 B psi2 92.333 D nx2 -0.920774 ny2 -0.264096 nz2 -0.287103 ! Script will include these End Molecule pyr1 at x1 y1 z1 rotated by psi1 about nx1 ny1 nz1 #include ./pyr.mom End Edit pyr1 #include ./pyr.axes End Polarizabilities for pyr1 ! Assumed that the pols are in the local-axes ! So they are read after axes are defined Read rank 1 #include ./pyr.pol End Limit rank 1 for +++ H1 H2 H3 H4 H5 N C1 C2 C3 C4 +++ C5 End Molecule pyr2 at x2 y2 z2 rotated by psi2 about nx2 ny2 nz2 #include ./pyr.mom End Edit pyr2 #include ./pyr.axes End Polarizabilities for pyr2 ! Assumed that the pols are in the local-axes ! So they are read after axes are defined Read rank 1 #include ./pyr.pol End Limit rank 1 for +++ H1 H2 H3 H4 H5 N C1 C2 C3 C4 +++ C5 End Pairs #include ./pyr2.pot End Units Bohr kJ/mol Switch Induce On Iterate On Options Induction Iterations 100 Convergence 1e-12 End Units Bohr kJ/mol Comment "Optimization" Optimize EVF ( This is the eigenvector following algorithm ) Search 0 ( search for minima ) Iterations 1000 ( number of iterations in algorithm ) Convergence 10e-3 ( in current units: kJ/mol ) Step 0.1 ( step in current units: Bohr ) Frequencies FILE pyr-n2 Animations Frames 50 Intensities Vectors End Comment "Energy calculation" Energy Details Analyse Finish
Most of this file is used just to define the PES, the molecular geometries, and the starting geometry of the pair of molecules. The functional part that deals with the geometry optimization is right at the end in the Optimize...End block. There we define the minimization algorithm to be used and various parameters associated with the algorithm. The Frequencies command enables calculation of vibrational frequencies of the system.
At the very end we perform an energy calculation and details of the energy are printed out. This is often useful.
Run the file using the command:
$ orient-4.8 < pyr2_geom_opt.ornt >& pyr2_geom_opt_model1.out &
The current version of Orient is 4.9.x but the data files this version of Orient expects are somewhat different in format from those we have here. In particular, the pyr.pol file needs to be in a different format. So we use Orient 4.8 in the above command. If you miss the -4.8 then you will get the more up-to-date version of Orient and it will not work on these data sets.
The output will be in file pyr2_geom_opt_model1.out. Have a look at this file. The output of the optimization section should be something like:
Optimization has converged after 4 iterations. Electrostatic energy: -12.881139 Repulsion energy: 15.231819 Induction energy: -2.709836 Dispersion energy: -15.754727 Total energy: -16.113883 kJ/mol Frequencies of the intermolecular modes (cm-1): 0.0073i 0.0069i 0.0000i 0.0000 0.0000 0.0005 15.9566 32.6517 43.1147 47.1893 54.6118 89.5815 Intensities written to file pyr-n2_intens.dat Molecules 1 and 2: intermolecular distance = 10.8705 bohr Orientations relative to 1->2 vector: Euler angles for molecule 1: 0.000 90.000 120.030 Euler angles for molecule 2: -180.000 90.000 -120.030
The code prints out the total energy of the system after optimization and the frequencies of the intermolecular models. For a minimum we should have all frequencies real. In this case we see that the imagiinary frequencies are all close to zero. That means we have a real minimum - though probably not he global minimum.