How to build isotactic PHB oligomers using mbuild.lib.recipes.polymer.Polymer

[janan94]

I am trying to build a 10-mer isotactic PHB oligomer using
mbuild.lib.recipes.polymer.Polymer.

I want the repeat unit to preserve stereochemistry, e.g. all repeat-unit
chiral centers have the same configuration.

example:

import mbuild as mb
from mbuild.lib.recipes.polymer import Polymer

monomer = mb.load("CC@HCC(O)=O", smiles=True)

I am unsure which atom indices should be used for the two ports

and whether Polymer preserves the stereochemistry from the SMILES.

chain = Polymer()
chain.add_monomer(
compound=monomer,
indices=[...],
separation=0.15,
replace=True
)
chain.build(n=10, sequence="A")

Questions:

1. Is this the recommended way to model PHB using mBuild?
2. Does Polymer preserve the stereochemistry from the chiral SMILES?
3. What is the correct way to define the two ports for ester linkage formation?
   Can anyone help me with this?

[chrisjonesBSU]

Hey @janan94. Sorry about the late response, it's been a busy week.

I'd suggest checking out the develop branch of mBuild to use new features like using smiles-string tagging (rather than finding particle indices) and using mbuild/path to create different chain configurations. Here's an example:

import mbuild as mb
from mbuild import Polymer
from mbuild.path import straight_line, Path

# Use SMILES-tagging with the {} brackets. You can put whatever you want here, no strict syntax
# In this case, the idea of tagging an atom in the SMILES string is to identify it in the polymer builder
monomer = mb.load("O{>}[C@H]CC{<}(=O)", smiles=True)
chain = Polymer()
# Pass in the tags you chose when loading the monomer from SMILES
chain.add_monomer(compound=monomer, head_tag=">", tail_tag="<", separation=0.15)
chain.build_from_path(path=straight_line(path=Path(), N=5, spacing=0.45))
# The spacing value is a rough estimate of the monomer-monomer center distance

These changes aren't in the conda-forge package, so to get the develop branch of mBuild there are a couple options:

1. Make a new mbuild environment:

git clone git@github.com:mosdef-hub/mbuild.git
cd mbuild
git checkout develop
conda env create -f environment.yml
conda activate mbuild
pip install .

Or if you already have an environment with mBuild installed you can:

git clone git@github.com:mosdef-hub/mbuild.git
cd mbuild
git checkout develop
conda activate your-environment
conda env update --file environment.yml --prune
pip install .

I'd recommend option 1 initially to have a new environment to test things out.

Let me or @CalCraven know if there are any issues.

[CalCraven]

Is this the recommended way to model PHB using mBuild?

See @chrisjonesBSU's response for other options. But the way you're doing it is certainly viable and should work if you can select the correct sites.

Does Polymer preserve the stereochemistry from the chiral SMILES?

Polymer will just place linearly the monomers, so the mb.load() method I think is where the chirality will be determined. This uses RDKit as the backend to get the compound coordinates, which I believe will read and configure based on the chirality. To confirm this, you can load two monomers with differing chiral centers, then use compound.visualize() to manually look at the atom of interest and verify they're both oriented differently and correctly. Make sure to do this in a jupyter notebook!

What is the correct way to define the two ports for ester linkage formation?

If you're using this older polymer builder, I recommend the old guess and check method. What I do, is something like this:

polymer_inds = [0,1] # randomly guess, but usually these are ordered by heavy atoms first as defined in the monomer smiles. 
sitesList = list(compound.particles())
for site in [sitesList[polymer_inds[0]], sitesList[polymer_inds[1]]]:
    site.name = "B" # set to be boron so it's easy to visualize
    site.element = ele.element_from_symbol("B")
compound.visualize()

Then look at the two particles, and manually use different values in polymer_inds to pick out which are the two hydrogens you'd want to reacted away.

You can replace the hydrogens if you're reacting away hydrogen sites, or you can use replace=False if you just want a bond to form connected to that atom.