AJMPublic/orient/optimization

Navigation:

• The Orient Program

Geometry optimization

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.

Data needed

• 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 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 output

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.