Tutorials

Tutorial 1: Force-field parameters from SMILES string

It is straightforward to obtain force-field parameters just from SMILES string by combining the functionality of metallicious with cgbind:

from cgbind import Linker, Cage
from metallicious import supramolecular_structure

# Create coordination file from cgbind
linker = Linker(smiles='C1(C#CC2=CC=CC(C#CC3=CC=CN=C3)=C2)=CC=CN=C1', arch_name='m2l4')
cage = Cage(linker, metal='Pd')
cage.print_xyz_file(filename="cage.xyz")

# Parametrize the cage, use AmberTools to parametrize organic linkers
cage = supramolecular_structure('cage.xyz', metal_charges={'Pd': 2}, LJ_type='uff')
cage.parametrize(out_coord='out.pdb', out_topol='out.top', prepare_initial_topology=True)

Tutorial 2: Parametrization of template for gallium cages (known issue and solution)

Metallicious uses educated guesses along the way. While we observe that this works in most cases, we also found one outlier: parametrization of the template for [Ga4L6]12- cages. To parametrize the template, metalliciuos uses the charge of the metal (user input) and guessed charge of the attached ligand fragment. The charge of the organic ligand is guessed by rdkit from the 3D structure of the atoms. While this works for most of the molecules, sometimes there might be ambiguity about the structure. In this example, the coordinating ligand can have two chemical structures:

metallicious incorrectly guesses the neutral charge of the ligand. However, this behaviour can be corrected by overwriting guessed charges:

from metallicious import supramolecular_structure
cage = supramolecular_structure('nonbonded.pdb', {'Ga': (3, 1)},
                                topol = 'noncovalent_complex.top',
                                LJ_type='uff', search_library=False)
cage.extract_unique_metal_sites()
cage.unique_sites[0].ligand_charges = [-2, -2, -2] #this specifies ligands' charges
cage.parametrize()

While it is worth noting this, this template is part of the default library, so it is already parametrized.

Tutorial 3: Using custom template

One might decide that they parametrized the template using a specific methodology (explicit solvent, or diffrent than usuall level of theory) and did not save it to the library. One might use this template on purpose for specific application

cage = supramolecular_structure(f'cage.xyz',
                                metal_charges={'Pd': 2},
                                LJ_type = 'uff',
                                search_library=False)
cage.extract_unique_metal_sites() # extracts template structures (needed to detect sites)
cage.sites[0].fp_coord_file = 'selected_template.pdb' # loads coordinates of the template
cage.sites[0].fp_topol_file = 'selected_template.top' # loads force-field parameters of the template
cage.sites[0].load_fingerprint() # loads template into the class
cage.sites[0].set_cutoff() # reads the radius of the template
cage.unique_sites = [] # overwrites list with templates to parameterize (no template parametrization needed)
cage.parametrize(out_coord='saturated_template.pdb', out_topol='saturated_template.top')

Tutorial 4: Custom template library

As default metallicious has a library of the templates parametrized using D3BJ-PBE0/def2-SVP and they are saved in directory of metallicious in library folder. One might however opt for higher level of theory or include implicit solvent effect. The easiest is to create new directory for library of templates. This requires only change of the library_directory, as there are no templates inside, we need to parametrize them using new method:

from metallicious import supramolecular_structure
cage = supramolecular_structure('nonbonded.pdb', {'Pd': (2, 1)},
                                topol = 'noncovalent_complex.top',
                                LJ_type='uff'
                                library_path=f'/path/to/new/library/',
                                keywords = ['CPCM(Water)', 'PBE0', 'D3BJ', 'def2-SVP', 'tightOPT', 'freq'])
cage.parametrize()

Tutorial 5: GROMACS tutorial for MD simulations of metallorganic structures from scratch

This tutorial will guide you in performing an MD simulation from scratch using GROMACS. We will simulate a Fujita’s Pd4L6 cage in water. The files needed for this tutorial are available here.

Firstly, create the topology (force-field parameters) files from .xyz file using metallicious:

from metallicious import supramolecular_structure
cage = supramolecular_structure('cage_start.pdb', metal_charges={'Pd': 2}, LJ_type='uff')
cage.parametrize(out_coord='0_cage.pdb', out_topol='topol.top', prepare_initial_topology=True)

Execute the script:

python script.py

With the topology ready, we will set up the system. Modify the size of the simulation box to accommodate the cage (we use 2 nm distance from the cage to the wall):

gmx editconf -f 0_cage.pdb -d 2 -o 1_box.pdb

Solvate the system with water (here, we use a 3-site water model):

gmx solvate -cp 1_box.pdb -cs spc216.gro -o 2_solv.gro -p topol.top

Change the topology files. Open the “topol.top” file and change the header:

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             0.5          0.833333

to force-field parameters, which include atom types of water:

#include "amber99sb-ildn.ff/forcefield.itp"

This will include all LJ parameters for atoms (including solvent and ions). Next, add information about the force-field parameters for water and ions just under the section “[ atomtypes ]” (before the first [ moleculetype ]). In this case, we will use SPC/E water model:

#include "amber99sb-ildn.ff/spce.itp"
#include "amber99sb-ildn.ff/ions.itp"

Now, we will add counter ions. Firstly, create a .tpr file:

gmx grompp -f em.mdp -c 2_solv.gro -p topol.top -o ions.tpr -maxwarn  2

Replace the water molecules with ions:

gmx genion -s ions.tpr -o 3_ions.gro -p topol.top -pname NA -nname CL -neutral

Select “SOL” to replace some of bulk water molecules with the counter ions to neutralize the cage charge. Minimize the created system (to “remove bad contacts”) by making a .tpr file and performing minimization.

gmx grompp -f em.mdp -c 3_ions.gro -p topol.top -o 4_em.tpr
gmx mdrun -v -deffnm 4_em

Equilibrate system using NVT and NPT ensemble:

gmx grompp -f nvt.mdp -c 4_em.gro -p topol.top -o 5_nvt.tpr
gmx mdrun -v -deffnm 5_nvt

gmx grompp -f npt.mdp -c 5_nvt.gro -p topol.top -o 6_npt.tpr
gmx mdrun -v -deffnm 6_npt

Finally, perform production run (ideally done on a high performance cluster (HPC) with GPUs):

gmx grompp -f run.mdp -c 6_npt.gro -p topol.top -o 7_run.tpr
gmx mdrun -v -deffnm 7_run

As a result, you should obtain the final frame 7_run.gro and trajectory 7_run.xtc, which you can visualise using a molecular visualization program (e.g., VMD). For your convenience, here are the final files compliled.

Tutorial 6: Host-guest complex parametrization

Parametrization of the host-guest complex involves separate parametrization of the host and guest. In this example, we will parametrize benzoquinone bound in the [Pd2L4]4+ cage:

Firstly, extract the host and guest structures into separate coordination files. While visualization programs like Avogadro or PyMOL can be used for this task, a more straightforward approach is to use OpenBabel:

obabel -ixyz start.xyz -oxyz -O separate.xyz -m --separate

This command will create two separate files: separate1.xyz and separate2.xyz. In this case, separate1.xyz corresponds to the guest, and separate2.xyz to the cage.

The guest molecule can be parametrized using the antechamber interface in the metallicious package:

from metallicious.antechamber_interface import antechamber
antechamber("separate1.xyz", "ligand.top", charge=0)

This will produce two files: ligand.top (force-field parameters) and ligand.pdb (the coordination files).

Perform standard parametrization of the cage using metallicious:

from metallicious import supramolecular_structure
cage = supramolecular_structure("separate2.xyz", metal_charges={'Pd':2 }, LJ_type='merz-opc')
cage.prepare_initial_topology()
cage.parametrize(out_coord=f'cage_out.pdb', out_topol=f'cage_out.top')

This will generate the coordination file (cage_out.pdb) and force-field files (cage_out.top) for the cage. (Remember that metallicious reorders atoms, so for the next steps, you need to use cage_out.pdb, not separate2.xyz).

To merge the coordinate files for the guest and cage, use the MDAnalysis package:

import MDAnalysis
syst1 = MDAnalysis.Universe("ligand.pdb")
syst2 = MDAnalysis.Universe("cage_out.pdb")
syst_sum = MDAnalysis.Merge(syst1.atoms, syst2.atoms)
syst_sum.atoms.write('complex_out.pdb')

The output file, complex_out.pdb, will contain the combined coordinates of the host-guest complex.

To merge the force-field parameters, use parmed:

import parmed as pmd
topol1 = pmd.load_file("ligand.top")
topol2 = pmd.load_file("cage_out.top")
topol_sum = topol1 + topol2
topol_sum.write("complex_out.top")

This will produce complex_out.top file containing the combined parameters of the guest and cage.

At the end of this process, you should obtain complex_out.pdb (coordination file) and complex_out.top (force-field parameters), which can be used for molecular dynamics with GROMACS. You can save them in different formats if you want to use them with a different molecular dynamics engine (for example, force-field with *.prmtop and coordination file with *.inpcrd, which can be used with AMBER).