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

  1. Filippini and A. Gavezzotti, Acta Cryst. B49 (1993) 868-880

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:

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