## page was renamed from AJMGrpOnly/orient/basin-hopping <> Navigation: * [[AJMPublic/orient|The Orient Program]] * [[AJMPublic/camcasp|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 minimisation 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 [[ajm/potentials/pyridine-model1|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. {{{#!wiki warning Warning The data for Model(1) can be read by Orient 4.8 and earlier versions only. }}} ==== Pyridine dimer : Orient 4.8 commands ==== {{{#!wiki warning Warning This command file is suitable for Orient 4.8. It will not work for Orient 4.9. }}} * [[attachment: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: * [[attachment: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 }}} 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 * [[attachment:partial output]] {{{ 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.