Creating a Solvation Box

In this tutorial we will cover:

• how we can use the Extension’s PDBFixer integration to solvate a molecule

• how we may use BuildAMol to create a solvation box from scratch

When performing molecular dynamics (MD) the question of solvent can be addressed in two ways. Either there are explicit solvent molecules added to the simulation system or the solvent effects are implicitely included when computing the energy heuristic.

There are great software tools around that can generate solvation boxes around molecules of interest, often working in synergy with the MD engine you are going to use later on. As such, it is always a good idea to use a dedicated tool that can do exactly what you need. BuildAMol does not itself concern itself with MD or solvation processes, and therefore does not offer a built-in method to create solvation boxes. However, there are still ways we can construct such systems with a little creativity. The purpose of this tutorial is two-fold. For once, we will see how we can create a bounding box using an established software-port (PDBFixer), and secondly, we will see a from-scratch construction of a (admittedly somewhat crystalline) bounding box that may not be very useful in practice but hopefully illustrates the process of creating multi-molecule systems. With that said, let’s start!

import plotly
plotly.offline.init_notebook_mode()

Let’s start by quickly making some molecule. For instance a glycan…

import buildamol as bam

bam.load_sugars()

glc = bam.get_compound("GLC")
nag = bam.get_compound("NAG")

# assemble some glycan (using pre-installed glycosydic linkages)
mol = (glc % "14bb" * 3) @ 2 % "12bb" + nag
mol.optimize()

mol.show()

/opt/anaconda3/lib/python3.12/site-packages/jupyter_client/session.py:203: DeprecationWarning:

datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).

Using PDBFixer

One of the most convenient tools to use when it comes to solvation boxes is the Python package PDBFixer. As of version 1.2.9 BuildAMol includes an extension to use PDBFixer to construct watery solvation boxes including ions around molecules.

The usage of the extension is quite straightforward:

1. import the necessary module

2. run the solvate function

3. that’s it…

# at the moment the only supported backend is PDBFixer, so we can directly import the (standardized) "solvate" function
from buildamol.extensions.molecular_dynamics.solvate import solvate

# let's solvate the glycan in a water box with a 10 Å buffer
solvated = solvate(mol, box_padding=10, ionic_strength=5)

# now we can visualize the solvated system
solvated.show()

/opt/anaconda3/lib/python3.12/site-packages/jupyter_client/session.py:203: DeprecationWarning:

datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).

And there we have our solvated molecule! Notice how the solvation process added two new chains (B=water, C=ions). Also notice that none of the newly added water residues have bonds (which is irrelevant if you just want to do MD afterward but keep it in mind)! At this point solvated is a new Molecule instance that can be further modified in any way you like.

Building a Solvation Box from scratch

Now the “interesting” part of the tutorial begins ;-) As mentioned before, BuildAMol does not offer any inherent functionality to construct solvation boxes so we need to think of our own setup to do so. Essentially we will be using a combination of a for loop, as well as the move_to and merge methods to construct a multi-molecule system.

For starters, let’s demonstrate how we can generate a system of multiple water molecules arranged in a static grid like so:

# get a water molecule
water = bam.molecule("HOH")

# make an empty molecule which will be the system of water molecules
system = bam.Molecule.empty()

# now iterate over a grid of points and add a water molecule at each point
stride = 3
size = 25
for x in range(0, size, stride):
    for y in range(0, size, stride):
        for z in range(0, size, stride):

            # make a new water molecule
            _water = water.copy()
            # move it to the current position
            _water.move_to([x, y, z])
            # add it to the system
            system.merge(_water)

# to make everything a bit tidier, we will squash the system down to a single chain
# (if we don't do this every water molecule will have it's own chain...)
system.squash_chains()

# show the system
system.show()

/opt/anaconda3/lib/python3.12/site-packages/jupyter_client/session.py:203: DeprecationWarning:

datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).

So that’s how we can make simple solvation boxes with BuildAMol!

Now we can start thinking of how we might create a solvation box around our glycan molecule. We would like to avoid placing water molecules on top of the actual glycan structure so let’s define a “no go zone” based on the glycan atoms’ position…

import numpy as np

# move the glycan to the center of the system
# (the two are still completely separate but the glycan is in the correct spacial position now)
mol.move_to((size / 2, size / 2, size / 2))

# get the atom coordinates
coords = np.array([i.coord for i in mol.get_atoms()])

# now we define a function that will check if a given point is too close to a glycan atom
def is_too_close(coord, threshold=3.5):
    # calculate the distance to each atom
    distances = np.linalg.norm(coords - coord, axis=1)
    # check if any of the distances are less than 3.5 angstroms
    return np.any(distances < threshold)

# now we can iterate over the water molecules and remove those that are too close to the glycan
for water in system.get_residues():

    # check if the water molecule is inside the glycan
    if is_too_close(water.coord):
        # if it is, remove it
        system.remove_residues(water)

Now before we add the glycan to the system, let’s quickly visualise everything, to see if this all worked as expected.

# and now we can visualise the system
v = bam.utils.visual.MoleculeViewer3D()
v.draw_edges(*mol.get_bonds(), linewidth=4, showlegend=False)

# draw each water molecule as a small blue point
for residue in system.get_residues():
    v.draw_point(residue.serial_number, residue.coord, color="blue", size=5, opacity=0.4, showlegend=False)

v.show()

Now admittedly, the water molecules in the current box are very crystalline. There are a number of things we could do to change that to a achieve a more “life like” setup. For example, PDBFixer uses OpenMM to get a better placement and orientation of the water molecules. This does require a quick MD run, however. The most low-tech option would be a random orientation of each water molecule in the system. To do so we can add a rotate call to our code when generating the solvation box like so (notice a number of small changes compared to the code above):

water = bam.molecule("HOH")
system = bam.Molecule.empty()

stride = 2
size = 10
center = mol.center_of_geometry
for x in range(-size, size, stride):
    for y in range(-size, size, stride):
        for z in range(-size, size, stride):

            _x = x + center[0]
            _y = y + center[1]
            _z = z + center[2]

            # make a new water molecule
            _water = water.copy()
            # move it to the current position
            _water.move_to([_x, _y, _z])

            # check if the water molecule is inside the glycan
            if is_too_close(_water.center_of_geometry, threshold=2.5):
                continue

            # rotate it to a random orientation
            _water.rotate(angle=np.random.rand()*360, axis=np.random.rand(3))

            # add it to the system
            system.merge(_water)

system.squash_chains()

# merge the glycan and the water system
system.merge(mol)
system.show()

/opt/anaconda3/lib/python3.12/site-packages/jupyter_client/session.py:203: DeprecationWarning:

datetime.datetime.utcnow() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.now(datetime.UTC).

Alright, and there we have a more organic-looking waterbox. Of course, no ions are present yet, but at this point you probably have an inkling of how to add these. Assuming we are happy with our results so far, the only thing left to do is save it to a file for use with our favourite MD engine…

# now we can save the system to a file
system.to_pdb("./files/glycan_with_waterbox.pdb")

This is it for this tutorial. Hopefully you could see that, using a little creativity, we can use BuildAMol even for tasks that the library is not explicitly designed to do. Thanks for checking out this tutorial and good luck with your research!