Generate a Pore

This section illustrates how to create a cylindric pore and functionalize the surface as described in the supporting information of Kraus et al. (2020). The PoreMS package can be imported using

import porems as pms

An empty pore will be generated by simply initializing a new object with the desired properties.

pore = pms.PoreCylinder([8, 8, 10], 4.8, 5.5)

Once the initialization is done, it is possible to attach molecules to the surface inside and outside of the pore. For this, attachment methods are implemented. In the following example the pore system used in Ziegler et al. (2019) will be replicated. Molecule structures and topologies are available on the corresponding repository.

The surface functionalization is divided in two parts, defining the molecules, and attaching them to the surface. In the present case two catalyst molecules are needed on the inner pore surface whereas the remaining surface is filled with a high concentration of DMDMS and the outer surface is mostly functionalized with TMS.

First the molecules need to be defined in the needed format. For this, a molecule class is provided which can be visualized as a molecule construction kit. Multiple possibilities exist to define the molecules. As an example, three methods will be illustrated.

Trimethylsilyl or for short TMS is provided in the PoreMS package as a generic molecule and can be simply imported

tms = pms.gen.tms()

The catalyst molecule can be constructed by using a structural file as an input. Herefore, the structure file catabm.gro provided in the mentioned repository will be used.

catalyst = pms.Molecule("catalyst", "CAT", "catabm.gro")

Finally, DMDMS will be constructed as a new molecule.

# Initialize new molecule
dmdms = pms.Molecule("dmdms", "DMS")

# Define dictionaries for bond length and angles
b = {"sio": 0.155, "sic": 0.186, "ch": 0.109 ,"co": 0.123}

# Build molecules
dmdms.add("Si", [0, 0, 0])
dmdms.add("O", 0, r=b["sio"])
dmdms.add("Si", 1, r=b["sio"])

for i in range(3):
    if i==0:
        dmdms.add("O", 2, r=b["sio"], theta=50, phi=120*i)
    else:
        dmdms.add("C", 2, r=b["sic"], theta=50, phi=120*i)

dmdms.add("C", 3, r=b["co"])

# Add hydrogens
for i in range(4, 6+1):
    for j in range(3):
        dmdms.add("H", i, r=b["ch"], theta=30, phi=120*j)

Once the molecules are defined, they can be attached on the surface. Since the functionalization is done iteratively, special placement, and molecules with a lower concentration on the surface have a higher priority and should be attached first.

Since the catalyst molecules are placed far enough apart so they do not interact or influence each other. Therefore, a point symmetrical placement is chosen, for which a special attachment function is used.

pore.attach_special(mol=catalyst, mount=37, axis=[34, 22], amount=2, symmetry="point")

The other two molecules will be attached using the conventional attachment function

pore.attach(mol=dmdms, mount=0, axis=[1, 2], amount=0.3, site_type="in", inp="molar")
pore.attach(mol=tms, mount=0, axis=[1, 2], amount=0.3, site_type="ex", inp="molar")

After finishing the surface functionalization, the pore needs to be finalized, which fills empty binding sites with silanol groups creating the final structure

pore.finalize()

In order to show the properties of the generated pore, use the table function

print(pore.table())

This returns a pandas data frame of pore properties and allocation.

At this point the pore generation is completed and what is left is converting the programs data structure into a readable file-format using the functionalities of the Store class. For this a store function is provided that creates a structure file in the GROMACS format, a main topology containing the number of atoms, a topology for the basic surface groups and grid atoms and a pickle file of the pore object

pore.store()

To sum it up, the complete code is as follows

import porems as pms


# Create TMS molecule
tms = pms.gen.tms()

# Create catalyst molecule
catalyst = pms.Molecule("catalyst", "CAT", "catabm.gro")

# Create DMDMS molecule
## Initialize new molecule
dmdms = pms.Molecule("dmdms", "DMS")

## Define dictionaries for bond length and angles
b = {"sio": 0.155, "sic": 0.186, "ch": 0.109 ,"co": 0.123}

## Build molecules
dmdms.add("Si", [0, 0, 0])
dmdms.add("O", 0, r=b["sio"])
dmdms.add("Si", 1, r=b["sio"])

for i in range(3):
    if i==0:
        dmdms.add("O", 2, r=b["sio"], theta=50, phi=120*i)
    else:
        dmdms.add("C", 2, r=b["sic"], theta=50, phi=120*i)

dmdms.add("C", 3, r=b["co"])

## Add hydrogens
for i in range(4, 6+1):
    for j in range(3):
        dmdms.add("H", i, r=b["ch"], theta=30, phi=120*j)


# Initialize pore
pore = pms.PoreCylinder([8, 8, 10], 4.8, 5.5)

# Attach Catalyst
pore.attach_special(mol=catalyst, mount=37, axis=[34, 22], amount=2, symmetry="point")

# Attach surface molecules
pore.attach(mol=dmdms, mount=0, axis=[1, 2], amount=0.3, site_type="in", inp="molar")
pore.attach(mol=tms, mount=0, axis=[1, 2], amount=0.3, site_type="ex", inp="molar")

# Finalize pore
pore.finalize()

# Store pore
pore.store()


# Print pore properties
print(pore.table())