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:

Data needed

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



      Sites     100 polarizable     100
      S-functions 50000
      Alphas 50000
      Parameter-sets 50000
      Pairs 100000
      Molecules 10

      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

! These variables define the starting geometry for the pyridine dimer
! You may alter this.

  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

Molecule  pyr1 at  x1  y1  z1  rotated by psi1  about nx1  ny1  nz1
  #include ./
Edit pyr1
   #include ./pyr.axes
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
   Limit rank 1 for +++
      H1  H2  H3  H4  H5  N   C1  C2  C3  C4   +++

Molecule  pyr2 at  x2  y2  z2  rotated by psi2  about nx2  ny2  nz2
  #include ./
Edit pyr2
   #include ./pyr.axes
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
   Limit rank 1 for +++
      H1  H2  H3  H4  H5  N   C1  C2  C3  C4   +++

   #include ./pyr2.pot

Units Bohr kJ/mol

Switch Induce On Iterate On
  Induction Iterations 100  Convergence 1e-12

Units Bohr kJ/mol

Comment "Optimization"
  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

Comment "Energy calculation"
Energy Details


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

AJMPublic/orient/optimization (last edited 2021-04-14 12:56:43 by apw109)