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

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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

sas_211 = SlabAdsorptionSites(slab_211, surface='fcc211')

# 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()
model = maxdcp(atoms, adsorbate_species='CO', coverage=0.67, 
               adsorption_sites=sas_211)
view(model) 

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

2. 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()
model = mindcp(atoms, adsorbate_species='CO',
               min_adsorbate_distance=1.8,
               adsorption_sites=sas_211)
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()
    model = mindcp(atoms, adsorbate_species='CO', 
                   min_adsorbate_distance=1.8,
                   adsorption_sites=sas_211)
    images.append(model)
view(images) 

For more usage of the ACAT package, I strongly recommend this Jupyter notebook tutorial.