lukasturcani / stk

A Python library which allows construction and manipulation of complex molecules, as well as automatic molecular design and the creation of molecular databases.
https://stk.readthedocs.io/
MIT License
252 stars 47 forks source link

Connecting metal-organic cages axial position #386

Closed Jaloca2 closed 3 years ago

Jaloca2 commented 3 years ago

First of all, thank you for the excellent tutorials and stk code. I'm a beginner at that, interested in material science mainly in extended frameworks. I can form metal-organic cages based on paddle-wheel units (M2L4, M3L6). However, I cannot extend the cage in a coordination network using the axial position (my purpose is to use the cage as bb). I tried to put a Br\NH2 functional group with materials studio, and after calling the .mol file and trying a reaction, I obtain an error (I think related to the metal valence). I tried also to construct a paddle-wheel unit using the range(6) such as octahedral metal complexes, but I obtain the same error at the same point. Do you know how I can form an extended network using metal-organic cages as bb through the axial positions? A 2D extended framework for example (using cof code)

Thank you in advance.

andrewtarzia commented 3 years ago

Hey @Jaloca2 - Welcome and thank you!

This is an interesting structure problem! I have previously written code that performs a similar construction, where a MOC is reacted at functional groups on its ligand backbone. I modified that code to do the reaction at the axial position (clearly, I have used a weird set of building blocks, but you can figure out how to replace them with what you want; let me know how it goes for you!).

Note that the way I have done this was to actually modify the Paddlewheel topology to have an axial position. i hope the code makes sense! I will suggest that the selection of COF topology may be a problem here, so let us know what you are trying to build and we can try to help!

@lukasturcani - what do you think of this solution??

import numpy as np
import stk
import stko

class AxialPaddlewheel(stk.metal_complex.MetalComplex):
    """
    Represents a metal complex topology graph.

    Two distinct metal building blocks, one with five functional
    groups, one with four.

    Binding ligand building blocks with two functional groups are
    required for this topology graph and axial ligand building
    blocks with one functional group are required.

    When using a :class:`dict` for initialization, a
    :class:`.BuildingBlock` needs to be assigned to each of the
    following numbers:
        | axial_metal: (0, )
        | other_metal: (1, )
        | ligands: (0, 1, 2, 3)
        | axial_ligand: (4, )

    See :class:`.MetalComplex` for more details and examples.

    """

    _metal_vertex_prototypes = (
        stk.metal_complex.vertices.MetalVertex(0, [0, 1, 0]),
        stk.metal_complex.vertices.MetalVertex(1, [0, -1, 0]),
    )
    _ligand_vertex_prototypes = (
        stk.metal_complex.vertices.BiDentateLigandVertex(2, [2, 0, 0]),
        stk.metal_complex.vertices.BiDentateLigandVertex(3, [0, 0, 2]),
        stk.metal_complex.vertices.BiDentateLigandVertex(4, [-2, 0, 0]),
        stk.metal_complex.vertices.BiDentateLigandVertex(5, [0, 0, -2]),
        stk.metal_complex.vertices.MonoDentateLigandVertex(6, [0, 2, 0]),
    )

    _edge_prototypes = (
        stk.Edge(
            id=0,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[0],
            position=[0.1, 0.5, 0],
        ),
        stk.Edge(
            id=1,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[0],
            position=[0.1, -0.5, 0],
        ),

        stk.Edge(
            id=2,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[1],
            position=[0, 0.5, 0.1],
        ),
        stk.Edge(
            id=3,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[1],
            position=[0, -0.5, 0.1],
        ),

        stk.Edge(
            id=4,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[2],
            position=[-0.1, 0.5, 0],
        ),
        stk.Edge(
            id=5,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[2],
            position=[-0.1, -0.5, 0],
        ),

        stk.Edge(
            id=6,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[3],
            position=[0, 0.5, -0.1],
        ),
        stk.Edge(
            id=7,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[3],
            position=[0, -0.5, -0.1],
        ),
        stk.Edge(
            id=8,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[4],
            position=[0, 1.5, 0],
        ),
    )

def main():
    rdkit_opt = stko.MetalOptimizer(
        metal_binder_distance=1.9,
        metal_binder_forceconstant=1.0e2,
        max_iterations=200,
    )
    fgs = [
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8]~[#1]',
            bonders=(1, ),
            deleters=(2, ),
        ),
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8X1]',
            bonders=(1, ),
            deleters=(),
        ),
    ]
    ligand_bb = stk.BuildingBlock.init_from_file('org_carb.mol', fgs)
    axial_bb = stk.BuildingBlock(
        smiles='C(N)O',
        functional_groups=(stk.AlcoholFactory(), ),
    )
    axial_bb.write('axial_bb.mol')
    print(ligand_bb, axial_bb)

    copper_atom = stk.BuildingBlock(
        smiles='[Cu+2]',
        functional_groups=(
            stk.SingleAtom(stk.Cu(0, charge=2))
            for i in range(5)
        ),
        position_matrix=np.array([[0, 0, 0]]),
    )

    # Use a custom Paddlewheel topology that has axial positions.
    copper_pw = stk.ConstructedMolecule(
        AxialPaddlewheel(
            metals={copper_atom: (0, 1)},
            ligands={
                ligand_bb: (0, 1, 2, 3),
                axial_bb: (4, ),
            },
            optimizer=stk.MCHammer(),
        )
    )
    copper_pw.write('copper_pw.mol')
    copper_pw = rdkit_opt.optimize(copper_pw)
    copper_pw.write('copper_pw_opt.mol')
    copper_pw = stk.BuildingBlock.init_from_molecule(
        copper_pw, (stk.BromoFactory(), )
    )
    ligand_bb = stk.BuildingBlock.init_from_file(
        'org.mol', (stk.BromoFactory(), )
    )
    print(copper_pw, ligand_bb)
    cage = stk.ConstructedMolecule(
        topology_graph=stk.cage.M2L4Lantern(
            building_blocks=(
                copper_pw,
                ligand_bb,
            ),
            optimizer=stk.MCHammer(num_steps=2000),
        )
    )
    cage.write('cage.mol')
    cage = rdkit_opt.optimize(cage)
    cage.write('cage_opt.mol')

    cage_bb = stk.BuildingBlock.init_from_molecule(
        molecule=cage,
        functional_groups=(stk.PrimaryAminoFactory(), ),
    )
    linker_bb = stk.BuildingBlock(
        smiles='C1(C(C(C1N)N)N)N',
        functional_groups=(stk.PrimaryAminoFactory(), ),
    )
    print(cage_bb, linker_bb)
    cof = stk.ConstructedMolecule(
        topology_graph=stk.cof.Square(
            building_blocks=(cage_bb, linker_bb),
            lattice_size=(2, 2, 1),
            # Setting scale_steps to False tends to lead to a
            # better structure.
            optimizer=stk.Collapser(scale_steps=False),
        ),
    )
    cof.write('cof.mol')

if __name__ == '__main__':
    main()
Jaloca2 commented 3 years ago

I use the following code to construct a M3L6 cage.

Cu_atom = stk.BuildingBlock(
  smiles='[Cu+2]',
  functional_groups=(
      stk.SingleAtom(stk.Cu(0, charge=2))
      for i in range(6)
  ),
  position_matrix=([0, 0, 0], ),
)

bi_1 = stk.BuildingBlock(
    smiles='C(=O)(O)Br',
    functional_groups=[
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8]~[#1]',
            bonders=(1, ),
            deleters=(2, ),
        ),
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8X1]',
            bonders=(1, ),
            deleters=(),
        ),
    ]
)
show_stk_mol(bi_1)

mcomplex = stk.ConstructedMolecule(
    stk.metal_complex.Paddlewheel(
        metals=Cu_atom,
        ligands=bi_1,
        reaction_factory=stk.DativeReactionFactory(
            stk.GenericReactionFactory(
                bond_orders={
                    frozenset({
                        stk.GenericFunctionalGroup,
                        stk.SingleAtom,
                    }): 9,
                },
            ),
        ),
        optimizer=stk.MCHammer(),
    )
)
show_stk_mol(bi_1)
show_stk_mol(mcomplex)

complex_bb = stk.BuildingBlock.init_from_molecule(mcomplex, (stk.BromoFactory(), ))

Hbdc = stk.BuildingBlock(
    smiles=(
        'C1=CC(=CC(=C1)Br)C#CC2=CC=C(C=C2)C#CC3=CC(=CC=C3)Br'
    ),
    functional_groups=[stk.BromoFactory()],
)

cageH = stk.ConstructedMolecule(
    stk.cage.M3L6(
        building_blocks=(
            complex_bb,
            Hbdc,
        ),
        optimizer=stk.MCHammer(),
    )
)
show_stk_mol(Hbdc)
show_stk_mol(cageH)
stk.MolWriter().write(cageH, 'cageM3L6.mol')

Then, my purpose is to use it as bb through the axial positions of one Rh atom (the external one) but doesn't work. To try to fix it, I have added a Br\NH2 in the Rh atoms (materials studio) and I have called it to react in a reaction (Honeycomb because the Cage has 3 positions to react and the ligand 2), but doesn't work again:

import stk

bb1 = stk.BuildingBlock.init_from_file('cageT_Br.mol', [stk.BromoFactory()])
bb2 = stk.BuildingBlock('C1=CC(=CC=C1Br)Br', [stk.BromoFactory()])

cof = stk.ConstructedMolecule(
    topology_graph=stk.cof.Honeycomb(
        building_blocks=(bb1, bb2),
        lattice_size=(3, 3, 1),
        # Setting scale_steps to False tends to lead to a
        # better structure.
        optimizer=stk.Collapser(scale_steps=False),
    ),
)

When I use the same code of cofs, but the functional groups are located in the linkers (for example in an M2L4 lantern) works perfectly.

import stk

bb1 = stk.BuildingBlock.init_from_file('cageCHO.mol', [stk.AldehydeFactory()])
bb2 =  stk.BuildingBlock('C1=CC(=CC=C1N)N', [stk.PrimaryAminoFactory()])

cof = stk.ConstructedMolecule(
     topology_graph=stk.cof.PeriodicKagome(
        building_blocks=(bb1, bb2),
        lattice_size=(3, 3, 1),
        # Setting scale_steps to False tends to lead to a
        # better structure.
        optimizer=stk.Collapser(scale_steps=False),

    ),
)
show_stk_mol(cof)
stk.MolWriter().write(cof, 'cof.mol')

Thank you so much for your help

andrewtarzia commented 3 years ago

This is difficult to read, please use 3` to generate nice code blocks. (the preview tab helps you see what your comment will look like)

I am glad it works though? As far as I can tell?

Jaloca2 commented 3 years ago

Sorry for the inconvenience. I have edited it with the 3'

andrewtarzia commented 3 years ago

Thank you! So you are getting the structure you were aiming for? If so, @lukasturcani can close this issue.

Please feel free to ask more questions or join our discord if you want to discuss stk usage!

Jaloca2 commented 3 years ago

No, I cannot obtain the desired structure. Because what I want is to connect the M3L6 cage through the axial positions (to achieve a 2D extended network with a ditopic linker).

andrewtarzia commented 3 years ago

The first issue with your code is that the Kagome COF requires a 4 functional group BB and a 2 functional group BB but the M3L6 with FGs on the axial positions of the metals will be a 3 functional group BB.

Below, I use the Honeycomb Cof topology with the M3L6 cage. Maybe copy that into the jupyter notebook and build from there?

import numpy as np
import stk
import stko

class AxialPaddlewheel(stk.metal_complex.MetalComplex):
    """
    Represents a metal complex topology graph.

    Two distinct metal building blocks, one with five functional
    groups, one with four.

    Binding ligand building blocks with two functional groups are
    required for this topology graph and axial ligand building
    blocks with one functional group are required.

    When using a :class:`dict` for initialization, a
    :class:`.BuildingBlock` needs to be assigned to each of the
    following numbers:
        | axial_metal: (0, )
        | other_metal: (1, )
        | ligands: (0, 1, 2, 3)
        | axial_ligand: (4, )

    See :class:`.MetalComplex` for more details and examples.

    """

    _metal_vertex_prototypes = (
        stk.metal_complex.vertices.MetalVertex(0, [0, 1, 0]),
        stk.metal_complex.vertices.MetalVertex(1, [0, -1, 0]),
    )
    _ligand_vertex_prototypes = (
        stk.metal_complex.vertices.BiDentateLigandVertex(2, [2, 0, 0]),
        stk.metal_complex.vertices.BiDentateLigandVertex(3, [0, 0, 2]),
        stk.metal_complex.vertices.BiDentateLigandVertex(4, [-2, 0, 0]),
        stk.metal_complex.vertices.BiDentateLigandVertex(5, [0, 0, -2]),
        stk.metal_complex.vertices.MonoDentateLigandVertex(6, [0, 2, 0]),
    )

    _edge_prototypes = (
        stk.Edge(
            id=0,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[0],
            position=[0.1, 0.5, 0],
        ),
        stk.Edge(
            id=1,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[0],
            position=[0.1, -0.5, 0],
        ),

        stk.Edge(
            id=2,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[1],
            position=[0, 0.5, 0.1],
        ),
        stk.Edge(
            id=3,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[1],
            position=[0, -0.5, 0.1],
        ),

        stk.Edge(
            id=4,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[2],
            position=[-0.1, 0.5, 0],
        ),
        stk.Edge(
            id=5,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[2],
            position=[-0.1, -0.5, 0],
        ),

        stk.Edge(
            id=6,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[3],
            position=[0, 0.5, -0.1],
        ),
        stk.Edge(
            id=7,
            vertex1=_metal_vertex_prototypes[1],
            vertex2=_ligand_vertex_prototypes[3],
            position=[0, -0.5, -0.1],
        ),
        stk.Edge(
            id=8,
            vertex1=_metal_vertex_prototypes[0],
            vertex2=_ligand_vertex_prototypes[4],
            position=[0, 1.5, 0],
        ),
    )

def main():
    rdkit_opt = stko.MetalOptimizer(
        metal_binder_distance=1.9,
        metal_binder_forceconstant=1.0e2,
        max_iterations=200,
    )
    fgs = [
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8]~[#1]',
            bonders=(1, ),
            deleters=(2, ),
        ),
        stk.SmartsFunctionalGroupFactory(
            smarts='[#6]~[#8X1]',
            bonders=(1, ),
            deleters=(),
        ),
    ]
    ligand_bb = stk.BuildingBlock.init_from_file('org_carb.mol', fgs)
    axial_bb = stk.BuildingBlock(
        smiles='C(N)O',
        functional_groups=(stk.AlcoholFactory(), ),
    )
    axial_bb.write('axial_bb.mol')
    print(ligand_bb, axial_bb)

    copper_atom = stk.BuildingBlock(
        smiles='[Cu+2]',
        functional_groups=(
            stk.SingleAtom(stk.Cu(0, charge=2))
            for i in range(5)
        ),
        position_matrix=np.array([[0, 0, 0]]),
    )

    # Use a custom Paddlewheel topology that has axial positions.
    copper_pw = stk.ConstructedMolecule(
        AxialPaddlewheel(
            metals={copper_atom: (0, 1)},
            ligands={
                ligand_bb: (0, 1, 2, 3),
                axial_bb: (4, ),
            },
            optimizer=stk.MCHammer(),
        )
    )
    copper_pw.write('copper_pw.mol')
    copper_pw = rdkit_opt.optimize(copper_pw)
    copper_pw.write('copper_pw_opt.mol')
    copper_pw = stk.BuildingBlock.init_from_molecule(
        copper_pw, (stk.BromoFactory(), )
    )
    ligand_bb = stk.BuildingBlock.init_from_file(
        'org.mol', (stk.BromoFactory(), )
    )
    print(copper_pw, ligand_bb)
    cage = stk.ConstructedMolecule(
        topology_graph=stk.cage.M3L6(
            building_blocks=(
                copper_pw,
                ligand_bb,
            ),
            optimizer=stk.MCHammer(num_steps=2000),
        )
    )
    cage.write('cage.mol')
    # cage = rdkit_opt.optimize(cage)
    # cage.write('cage_opt.mol')

    cage_bb = stk.BuildingBlock.init_from_molecule(
        molecule=cage,
        functional_groups=(stk.PrimaryAminoFactory(), ),
    )
    linker_bb = stk.BuildingBlock(
        smiles='C1=CC(=CC=C1N)N',
        functional_groups=(stk.PrimaryAminoFactory(), ),
    )
    print(cage_bb, linker_bb)
    cof = stk.ConstructedMolecule(
        topology_graph=stk.cof.Honeycomb(
            building_blocks=(cage_bb, linker_bb),
            lattice_size=(2, 2, 1),
            # Setting scale_steps to False tends to lead to a
            # better structure.
            optimizer=stk.Collapser(scale_steps=False),
        ),
    )
    cof.write('cof.mol')

if __name__ == '__main__':
    main()
andrewtarzia commented 3 years ago

Clearly, you will need to modify the building blocks to be the chemical structures you are targeting.

I hope this makes sense.

Jaloca2 commented 3 years ago

After several tries, it works!!! Amazing work. Thank you so much for your help! If I see you at a chemistry conference, I will invite you to a couple of beers, I promise.

andrewtarzia commented 3 years ago

Excellent! Always happy to help and I will get the second round (I will be virtually at EuroMOF if that is in your field?).

@lukasturcani please close this issue.