I want to model a reaction catalysed by periclase ($\ce{MgO}$) using DFT. I have a good guess on the transition state (TS) of the reaction that goes in gas phase/solvent (produced using MOPAC).

The proposal for this reaction goes by the Eley–Rideal mechanism and I have a good guess on how the reactant adsorbs.

Now comes the question: I would like to optimise the TS on the surface using DFT (either ORCA or NWChem), but in tractable time.

I am tempted to use a simple, fixed model for the surface, e.g., a set of fixed atoms in bulk arrangement. Can this work? How am I supposed to find a structure with just one imaginary frequency using a partially fixed model (maybe using some clever vibrational projections)?

Even though NEB could be applied, the implementation of it in NWChem does not work with fixed atomic coordinates. Besides, I am confident on the TS I found without the catalyst, so I would like to use that instead of starting all over again.

Edit: what about using a very small cluster model, say, 8-16 atoms, and let it relax? Do you think it could work (in the sense of giving a reliable TS in terms of geometry and energy)?

  • 2
    $\begingroup$ It will not be easy because the slab (that's what a cutout of a solid state structure is often called) would need to be optimized in order to supply a decent initial Hessian for the TS search. I doubt that any TS search algorithm deals well with fixed atoms. If you manage that, I would strongly suggest modelling the surrounding Mg/O by pseudopotentials and point charges. My old group did similar work for polarization properties: DOI:10.1002/cphc.201100521. $\endgroup$
    – TAR86
    Commented Jan 27, 2017 at 12:02
  • 1
    $\begingroup$ Several comments: mopac ts maybe qualitatively wrong; surface and vacuum ts-s are generally very different, therefore your approach can easily be wrong; orca cannot calculate slabs; slab calculations takes forever, but you cannot cut corners most cases, and so on $\endgroup$
    – Greg
    Commented Jan 27, 2017 at 16:21
  • 1
    $\begingroup$ With a cluster model, you may get an answer. As Greg points out: whether that has much to do with the reality, I have my doubts, but also lack experience: my old group mostly dealt with gas-phase TS. As a first step, you might check whether your cluster relaxes to a reasonable structure when surrounded by the artificial environment (point charges etc.), because without a relaxed geometry, I don't think you will get anywhere. Next, I would look at a trivial model reaction and test the remaining approach. Given those results, I would assess the feasibility. And don't forget dispersion/vdW. $\endgroup$
    – TAR86
    Commented Jan 27, 2017 at 17:04
  • 2
    $\begingroup$ I would question if there is a fast/easy/cheap solution that does not require experience (or even just a fast one) $\endgroup$
    – Greg
    Commented Jan 27, 2017 at 17:09
  • 2
    $\begingroup$ The standard scheme here is to take a slab and keep 2-4 atom layers outermost from the reaction center fixed in crystallographic positions while letting everything else to relax naturally. There is a plenty of works done in this fashion. The size of the model used (i.e. how much atomic layers counting from the reaction center) to be used is ideally a matter of trial and asymptotic analysis, but in my experience in published works usually 2-4 layers from reaction center are relaxed and 2-4 more are fixed. $\endgroup$
    – permeakra
    Commented Feb 20, 2017 at 12:43

1 Answer 1


To study surface reactions I recommend the growing string method (GSM) for surfaces. Here is a nice quote from the paper which developed it (Ref. 1)

GSM’s efficacy was confirmed by comparison with CI-NEB on an extensive set of reactions characteristic of modern surface chemistry studies. In these cases, GSM reduces the computational cost (in terms of gradient computations) by about 45% on average over CI-NEB. In addition to high efficiency, GSM has the advantage of operating in single-ended way to enable explorative study of chemical reactions.

Key to the success of GSM is the use of the efficient hybrid delocalized coordinate system which greatly improves optimization and interpolation steps in comparison to NEB which uses Cartesian coordinates.

The exact code to reproduce their work is not publicly available, although you may be able to request a copy. However you should be able to reproduce it and more, with this python version of the code that I have developed.

Also, as others have pointed out, you should consider a higher level of theory like periodic wavefunction DFT (e.g. as implemented in VASP) .


  1. Jafari, M, Zimmerman, P. M.. "Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method", J. Comput. Chem. 2017, 38, 645– 658. https://doi.org/10.1002/jcc.24720
  • $\begingroup$ Good answer and very interesting method. I'll take a look at it and accept your answer, thanks! $\endgroup$ Commented Apr 6, 2020 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.