# How to model a random distribution of molecules attached to a surface?

In the same aspect as some of my previous questions, I am interested in modeling molecules on different surfaces and to study the interfacial effect on them.

One of the most important problems that one can face in such aspect is how to create an initial guess of randomly oriented and distributed molecules over the material surface.

While it is relatively easy to calculate a single molecule on a surface and, if necessary, use Periodic Boundary Conditions to replicate such unit in the material plane, I am not aware in what is the best option to create a much denser 'packing' in a randomized distribution.

Any idea about software or methodology to set this initial structure?

• You want to use Packmol, definitely, which lets you create various random initial packings of molecules - whether on surfaces or in biomolecular simulations. One thing to remember is that molecules on surfaces are rarely perpendicular to the surface plane - there's usually a tilt angle. Commented May 11, 2020 at 12:56
• Thanks Geoff, I’ll take a look at such software. :) Commented May 11, 2020 at 12:57
• @GeoffHutchison, you should make that into an answer rather than a comment. Commented May 12, 2020 at 4:03

As mentioned by @GeoffHutchison, PACKMOL is probably the most used one for random dense packing of molecules.

Sometimes this might however not suit your needs, e.g. if the shape you want to fit molecules into is not supported by PACKMOL or you want to achieve a certain distance. I achieved good results using the surface tools of ASE. Since it's Python you can create random positions with restraints you choose. Of course you will lack the automatic overlap prevention PACKMOL offers.

• > Of course you will lack the automatic overlap prevention PACKMOL offers. i don't think this is true. see blmin = closest_distances_generator(atom_numbers=[Z], ratio_of_covalent_radii=0.5) passed to StartGenerator in this docs example Commented Sep 13 at 14:24

If your molecules form not too many types of adsorption complexes on your surfaces, then you can built a lattice model with SuSMoST ( http://SuSMoST.com/ ) and run a short Monte Carlo simulation at high temperature. High temperature Monte Carlo gives you essentially a random structure.

I might be late to the party, but I have recently developed a Python package ACAT that can help people with similar problems of building surface adsorption models. The package is interfaced with the popular atomistic modelling package ASE, and supports a wide variety of molecules, surfaces (e.g. fcc/bcc/hcp structures with various Miller indices) and even nanoparticles (fcc/icosahedron/decahedron).

Let me start by addressing your question specifically. You want to "model a random distribution of molecules attached to a surface", but oftentimes it is necessary to consider the following constraints to the randomness:

1. When molecules are attached to a surface, they prefer to only adsorb at the high-symmetry sites, i.e. ontop/bridge/hollow sites.
2. Even if one wants to generate a dense packing of molecules, it can become unphysical if the molecules are too close to each other (or even overlapping).

Therefore, to achieve what you want in a semi-automated fashion, you will need

1. a tool to generate the surface
2. a tool to automatically identify all high-symmetry adsorption sites on a given surface structure
3. a tool to add adsorbates randomly to the high-symmetry sites with certain constraints

A surface can be either generated by ASE or read in from a structure file, then we can use ACAT for the 2nd and 3rd needs. As an example, I will show here the Python code for generating a random coverage of CO molecules on an fcc(211) stepped Ni surface. First, we can generate the surface using ASE:

from ase.build import fcc211

slab_211 = fcc211('Ni', (6,3,4), vacuum=5.)


Then we can identify all high-symmetry adsorption sites on the fcc(211) surface using ACAT:

from acat.adsorption_sites import SlabAdsorptionSites

# You can also print out the information of each site
sites_211 = sas_211.get_sites()
for i, site in enumerate(sites_211):
print('Site {0}: {1}'.format(i, site))


Now that we want to generate random adsorbate coverage patterns, I have implemented 2 functions for this purpose:

1. max_dist_coverage_pattern uses a clustering algorithm to generate a coverage pattern that maximizes the minimum adsorbate-adsorbate distance given the positions of adsorption sites and the number of adsorbates to be added. For example, we can generate a random CO* coverage pattern on the Ni(211) surface with a coverage of 0.67 ML (monolayer) by:
from acat.build.adlayer import max_dist_coverage_pattern as maxdcp
from ase.visualize import view

atoms = slab_211.copy()
view(model)


And this is the top view of the final surface adsorption model that we obtain:

1. min_dist_coverage_pattern maximizes the density of the adsorbates when a minimum adsorbate-adsorbate distance constraint is given. Note that the number of adsorbates generated by this function is not fixed. For example, we can generate a random CO* coverage pattern on the Ni(211) surface with a minimum adsorbate-adsorbate distance of 1.8 Å:
from acat.build.adlayer import min_dist_coverage_pattern as mindcp
from ase.visualize import view

atoms = slab_211.copy()
view(model)


And this is the top view of the final surface adsorption model that we obtain:

If we want to generate multiple random coverage patterns, we can simply use a for loop:

images = []
for _ in range(100):
atoms = slab_211.copy()