Geometry Optimisation

Note

The scripts for this step of the tutorial can be found in mgo_surface/3_mgo_opt.

In this section we demonstrate an energy minimisation of the MgO cluster. The input script mgo_opt_energy.py is very similar in form to the single point calculation and you can refer back to the previous step for explanations of the identical parts.

To run a geometry optimisation we need to know which atoms are active, that is can be relaxed (corresponding to regions 1-3 of the ionic embedded cluster model). We can use the getRegion() method to get a list of active atoms from the input fragment in the same way as for the QM region.

# Get list of active atoms for optimisation
active_region = mgo.getRegion(1,2,3)
print("Active region is: ", active_region)

The QM/MM settings are the same as for the single point calculation with one exception: shl_maxcycles is reduced to 5. This makes sure that the geometry optimisation does not get stuck in a long sequence of shell relaxations at any particular step, while it should be sufficient to ensure that the shells are fully relaxed by the end of the geometry optimisation.

The energy minimisation is run by calling DL-FIND, the standard geometry optimisation module for ChemShell. An instance of the Opt() class is first set up containing DL-FIND options, and the optimisation run using the run() method.

# Run QM/MM optimisation
opt = Opt(theory=qmmm, 
          algorithm="lbfgs", trust_radius="energy",
          maxcycle=200, maxene=300, active=active_region)
opt.run()

Here the geometry optimisation has several variables:

  • algorithm is used to specify the optimisation algorithm used by DL-FIND. Here we use the default L-BFGS algorithm which is suitable for both small and large systems.

  • maxcycle determines the maximum number of optimisation cycles. The optimiser will cease after convergence or at the limit set here.

  • maxene determines the maximum number of energy and gradient calculations during the optimisation. The optimiser will cease if this is reached before convergence.

  • active specifies the atoms that will be relaxed during the optimisation. Here, we pass in the active_region list of region 1-3 atoms.

With NWChem running on 4 cores, the optimisation will run in approximately 45 minutes.

After the optimisation completes, the optimised structures are saved for later use in punch file format, and for visualisation in an XYZ file.

# Save optimised geometry
mgo.save('mgo_shells_opt.pun', 'pun')
mgo.save('mgo_shells_opt.xyz', 'xyz')

Finally, the result is validated against the pre-calculated optimised energy.