Diff for "AJMPublic/orient/basin-hopping"

Navigation:

• The Orient Program

• CamCASP

Basin-Hopping

The basin-hopping algorithm is a stochastic simulation technique to efficiently find local minima on a potential energy surface (PES). The PES for a complex containing small molecules, such as water, can be relatively simple and will normally support only a few minima. For such systems finding this minimum is not a problem and any minimization technique will work. However the PESs for larger molecules, or larger numbers of molecules can support large numbers of minima. For such systems more robust methods are needed to find these minima.

References

See the innumerable papers by David Wales for details of how this algorithm works. Here are a select few:

• Energy Landscapes, with applications to clusters, biomolecules and glasses, D. J. Wales, Cambridge (2003). This is a large book, but you will find many examples of the algorithm here.

• Potential energy and free energy landscapes, D. J. Wales and T. V. Bogdan, J. Phys. Chem. B, **110**, 20765-20776 (2006).

Sample Orient command file

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.

Pyridine dimer : Orient 4.8 commands

• pyridine_n2_basinhopping.ornt

pyridine_n2_basinhopping.ornt

Title pyr2 : Basin-hopping

UNITS BOHR

Parameters
      Sites     26 polarizable     26
      S-functions 100000
      Alphas 100000
      Parameter-sets 100000
      Pairs 100000
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

Molecule  pyridine1 at  0.0 0.0 0.0 rotated by 0.0 about 0.0 0.0 1.0
  #include ./Model1/pyr.mom
End
Edit pyridine1
  #include ./Model1/pyr.axes
End
Polarizabilities for pyridine1
   ! Assumed that the pols are in the local-axes
   ! So they are read after axes are defined
   Read rank 1
     #include ./Model1/pyr.pol
   End
   Limit rank 1 for +++
      H1  H2  H3  H4  H5  N   C1  C2  C3  C4   +++
      C5 
End

Molecule  pyridine2 at  0.0 0.0 10.0 rotated by 0.0 about 0.0 0.0 1.0
  #include ./Model1/pyr.mom
End
Edit pyridine2
  #include ./Model1/pyr.axes
End
Polarizabilities for pyridine2
   ! Assumed that the pols are in the local-axes
   ! So they are read after axes are defined
   Read rank 1
     #include ./Model1/pyr.pol
   End
   Limit rank 1 for +++
      H1  H2  H3  H4  H5  N   C1  C2  C3  C4   +++
      C5 
End

Units Bohr Hartree

Pairs
   #include ./Model1/pyr2.pot
End

Units Bohr kJ/mol

Switch Induce On Iterate On
Options
  Induction Iterations 60  Convergence 1e-12
End

Options
  bfgs
  Iterations 1000
  convergence 1d-5
  show none
  plot none
end

Time

Basin-hopping
  ! restart save.cfg
  ! Verbose
  Allocate
    minima 100
  end
  Temperature 500
  Seed 777  ( optional )
  Steps
    number 500
    !  A serious calculation would need many more steps than this
    maxdisp 0.25
    maxrot 30 ( degrees )
    blocks rotations 50 translations 50 both 100
  end
  Difference energy 1e-5 MI 0.5
  Container radius 7.0
  Quenching
    reject -100
  end
  Files
    xyz    pyr2_model1_n2_run1.xyz
    minima pyr2_model1_n2_run1.geom
  end
end

Time

Finish

Run this job using (I have used orient-4.8 as the executable. You need to make sure that you use the correct executable):

$ orient-4.8 < pyridine_n2_basinhopping.ornt > pyr_n2_BH.out & 

This can take a while as the above example uses 500 steps in the search process. The output is also quite long, but we can use the Cluster program that is part of the CamCASP suite to analyse it and provide a sensible summary of the results.

Analysing the basin-hopping results using Cluster

The Cluster program that is part of the CamCASP suite of codes can be used to perform an analysis of the results of the Basin-Hopping simulations. Here is a sample code:

• reorient-analyse.clt

reorient-analyse.clt

Title Re-orient and analyse clusters

Global
  Units Bohr Degrees kJ/mol
End

Orient
  Geom-File RDX_run1.geom
  Geom-File RDX_run2.geom  APPEND OFFSET 100
  Reorient  RDX2
  ! R-scalings 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.20 1.40 
  Sort
  Write  Format Orient All
  Analyse
    Sort
    Similarity  Moments  & Energy 
    Energy-Sigma 0.5 kJ/mol
    Moment-Sigma 40.0 AMU
    Similarity-cutoff 0.5
    Details
  End
End

Finish

In this example, the outputs of two Basin-Hopping simulations is analysed. Cluster uses the files containing the minima. These will be the files names with the command:

minima  <filename>

in the Basin-Hopping commands in the Orient command file. Here these are RDX_run1.geom and RDX_run2.geom. The APPEND command tells Cluster to append the contents of the second file after those of the first. The **OFFSET 100** command tells it to change the numbers of the minima in the second file by 100. This ensures that we have unique numbers to all minima. Of course, if your files contain more than 100 minima each then a larger offset will be needed.

Run this using

  $ cluster < reorient-analyse.clt > reorient-analyse-1.out

It will run in a few seconds or so. The output will be fairly large, but the important parts are at the bottom and should be something like

partial output

Configuration statistics
! Config List name : CfgList-1
! Configs read from file : RDX_run1.geom
UNITS  BOHR DEGREE KJ/MOL AMU BOHR^2
! ----------------------------------------------------------------------------------------------------------------
! INDEX          ENERGY         I_XX           I_YY           I_ZZ           NUM-MOLS       NUM-SIMILAR    SYMMETRY

! ----------------------------------------------------------------------------------------------------------------
         9     -38.673719    4816.657600   26521.382000   29029.103000              2              2    C1
        13     -32.596392    5097.135300   25781.361000   28326.087000              2              4    C1
         8     -31.943536    7273.543700   15751.536000   16833.182000              2              2    C1
        12     -31.369167    6174.222300   21108.671000   24421.995000              2              2    C1
         1     -31.160460    6116.310500   20960.708000   24786.811000              2              2    C1
         6     -30.967247    5696.724600   24537.264000   27691.237000              2              2    C1
         4     -30.870978    5711.535600   24567.950000   27727.494000              2              2    C1
         5     -28.023902    6687.209600   16258.763000   18142.056000              2              2    C1
        20     -25.506775    7418.221300   17500.675000   18939.431000              2              2    C1
        15     -25.460627    8009.331400   14364.374000   15967.696000              2              2    C1
         3     -25.366652    8256.744300   15029.776000   15387.661000              2              2    C1
        18     -25.236796    5352.739300   25368.417000   27793.813000              2              4    C1
        17     -25.146462    7353.113400   17708.961000   19113.039000              2              2    C1
        21     -24.510071    6882.945500   20990.439000   22027.654000              2              2    C1
        14     -24.469602    6896.643800   20927.123000   21948.175000              2              2    C1
         7     -23.014142    6294.826700   22029.994000   25021.199000              2              4    C1
         2     -22.727332    6003.548700   24192.692000   28466.521000              2              4    C1
        16     -22.052672    7586.187000   20946.833000   21977.863000              2              2    C1

End

Finish
 Exiting program cluster_operations

This listing is a summary of the minima that Cluster considers as unique. The fields include:

• INDEX : the unique index of the minimum
• ENERGY: energy in units specified
• I_XX, I_YY, I_ZZ: moments of inertia in principle axes.
• NUM-MOLS: Number of molecules in the cluster. Here only 2 in each as we searched for only dimers.
• NUM-SIMILAR: Number of minima that Cluster considers to be similar. We expect to find 3-5 similar structures for each of the low-energy minima.
• SYMMETRY: At present the code does not find the point-group symmetry, so it prints out 'C1' for all.

The clustering is done using the commands in the ANALYSE block:

  Analyse
    Sort
    Similarity  Moments  & Energy 
    Energy-Sigma 0.5 kJ/mol
    Moment-Sigma 40.0 AMU
    Similarity-cutoff 0.5
    Details
  End

Here the structure similarity is based on the moments of inertia and the energies. Each quantity is given a tolerance set by a standard deviation. This defines a similarity probability as a normal distribution with the specified width. Structures are deemed to be similar to a probability. The cutoff level is set by **Similarity-cutoff**. Here it is $0.5$ so structures with a similarity probability $p \ge 0.5$ are deemed to be similar.

The clustering will depend on the standard deviations chosen. This needs to be tuned to the problem of interest. It is quite likely that the **Energy-Sigma** is too large in the example above, and quite likely that the **Moment-Sigma** is too small. You need to experiment to decide.