Crystal Optimiser¶
The class ccdc.csp.crystal_optimiser.CrystalOptimiser
allows to perform geometry optimisations of crystal structures.
The methodology is that used in CSD-Theory’s Landscape Generator.
For minimisation of molecular geometry, see Descriptive molecule minimisation documentation.
The optimisers find and often converge in local minima in the potential energy surface. Users can execute additional calculations with different convergence tolerance values to further explore the potential energy surface and potentially identify more stable minima.
In this example, we will explore the various settings available to users and optimise a crystal structure.
Let us start by opening the crystal structure of entry ACSALA01.
>>> from ccdc import io
>>> csd_refcode = 'ACSALA01'
>>> csd_reader = io.EntryReader("CSD")
>>> crystal = csd_reader.entry(csd_refcode).crystal
We then instantiate the CrystalOptimiser class.
>>> from ccdc.csp.crystal_optimiser import CrystalOptimiser
>>> optimiser = CrystalOptimiser()
The optimisation settings can be accessed through the nested ccdc.csp.crystal_optimiser.CrystalOptimiser.Settings
class.
>>> optimiser_settings = optimiser.settings
If we wish to save the results of our optimisation, we need specify a working directory. A preferred file type (.mol2 or .cif) for the optimised crystal structure can also be specified. The name of the saved file will be assigned automatically as Optimised_[refcode].file_format. If a working directory is specified, a log file for the calculation will also be saved.
>>> optimiser.settings.working_dir = "home/jdoe/optimised_crystals"
>>> optimiser.settings.output_format = "mol2"
Next, we can select the desired force field to use for the optimisation.
>>> optimiser.settings.force_field = "CLP"
The full set of available force fields is:
Name |
Reference |
---|---|
Dreiding II |
S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys. Chem. 1990 (94) 8897 |
Gavezzotti |
|
Momany |
F.A. Momany, L.M. Carruthers, R.F. McGuire and H.A. Scheraga, J. Phys. Chem. 78 (1974) 1595 |
CLP |
Gavezzotti, A. (2011). New J. Chem. 35, 1360-1368 |
CSD-OPCS16 |
Cole, J.C., Groom, C.R., Read, M.G., Giangreco, I., McCabe, P., Reilly, A.M., Shields, G.P., Acta Cryst. (2016) B72, 530-541. |
We are almost ready to optimise our structure. We just need to specify a few options that determine the type of optimisation we want. For this example, we will perform a ‘Full’ optimisation, optimising molecular geometry as well as unit cell parameters:
>>> optimiser.settings.optimise_molecule = "Geometry"
>>> optimiser.settings.optimise_cell = True
- For the remaining settings, we can use the default values. These are:
ccdc.csp.crystal_optimiser.CrystalOptimiser.Settings.cell_limit
, the maximum % volume change allowed for the unit cell (10%)ccdc.csp.crystal_optimiser.CrystalOptimiser.Settings.convergence_tolerance
, the convergence limit for the calculation (0.8)ccdc.csp.crystal_optimiser.CrystalOptimiser.Settings.limiting_radius
, the limiting radius for calculating intermolecular interactions (30 Angstroms)ccdc.csp.crystal_optimiser.CrystalOptimiser.Settings.normalise_hydrogens
, whether hydrogen atoms positions are normalised to neutron diffraction data (True)
After our settings have been defined, we can run the optimiser on our crystal and return an instance of the nested class ccdc.csp.crystal_optimiser.CrystalOptimiser.Results
.
>>> optimisation_results = optimiser.optimise(crystal)
>>> optimisation_results.n_steps
17
We can access the scores for our optimised crystal structure, as well as :
>>> round(optimisation_results.starting_score, 2)
-117.72
>>> round(optimisation_results.final_score, 2)
-125.26
Note that the scores are negative, and indicate the stability of our crystal structure. This scores correspond to the lattice energy when both below conditions are met:
The Force field selected is “CLP”, “Dreiding II”, “Momany” or “Gavezzotti”
Optimise molecule in crystal is set to “Position”
We can compare the unit cell parameters of the optimised crystal structure with those of the starting one:
>>> optimisation_results.original_crystal.cell_lengths
CellLengths(a=11.43, b=6.591, c=11.395)
>>> optimisation_results.optimised_crystal.cell_lengths
CellLengths(a=11.345, b=6.540, c=11.259)
>>> optimisation_results.original_crystal.cell_angles
CellAngles(alpha=90.0, beta=95.68, gamma=90.0)
>>> optimisation_results.optimised_crystal.cell_angles
CellAngles(alpha=90.0, beta=99.699, gamma=90.0)
>>> round(optimisation_results.original_crystal.cell_volume, 2)
854.23
>>> round(optimisation_results.optimised_crystal.cell_volume, 2)
823.54
And finally, we can compare our crystal structures making use of the ccdc.crystal.PackingSimilarity
class.
>>> from ccdc.crystal import PackingSimilarity
>>> packing_sim = PackingSimilarity()
>>> packing_sim_results = packing_sim.compare(optimisation_results.original_crystal, optimisation_results.optimised_crystal)
>>> round(packing_sim_results.rmsd, 3)
0.249