|
⇤ ← Revision 1 as of 2021-02-04 12:29:32
Size: 10473
Comment:
|
Size: 10475
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 11: | Line 11: |
| =====References===== | ===== References ===== |
Contents
Navigation:
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.
warning
The data for Model(1) can be read by Orient 4.8 and earlier versions only.
Pyridine dimer : Orient 4.8 commands
warning This command file is suitable for Orient 4.8. It will not work for Orient 4.9.
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
FinishRun 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: <code | 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
- Sort
End
Finish
</code> 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 <code | 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
- 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
- 3 -25.366652 8256.744300 15029.776000 15387.661000 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
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
EndHere 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.