Coloring molecule with RESP charge on pymol #Pymol

To visualize 3D structure of molecule, PyMol is nice tool. What I would like to write on the post is how to visualize calculated RESP charge on pymol ;)
One idea is embed calculated RESP charge to b_factor of molecule pdb. PDB file can make easily from rdkit mol object.
A problem for me is how to edit b_factor of PDB file. I found good tool to solve it today, python package ‘BioPandas’!
http://rasbt.github.io/biopandas/ It can install from pypi or anaconda.
By using biopandas, user can handle pdb file like pandas dataframe.

At first I would like to show examples of biopandas. If you often working with PDB you feel it is very useful I think.
Following code is an example for loading PDB and retrieve protein and ligand data as pandas dataframe.

from biopandas.pdb import PandasPdb
pbdobj = PandasPdb('1atp.pdb')
#Check Atom data
pbdobj.df['ATOM'].head(3)
#Check ligand data
pbdobj.df['HETATM'].head(3)

Now I can access pdb file as pandas data frame, it means that it is easy to edit and analyze pdb data! More examples are described in official site. So I move to next.

Following example is …
Calculate RESP charge with Psikit and then make pdf file from rdkit mol object and finally replace the b_factor to calculated RESP charge.
At frist load packages.

import os
import sys
from biopandas.pdb import PandasPdb
import pandas as pd
sys.path.append('path for /psikit')
from psikit import Psikit

Then define the function which make pdb file.

def make_pdb_with_Resp(smi, fname=None):
    pk = Psikit()
    pk.read_from_smiles(smi)
    pk.optimize()
    charge = pk.resp_charge
    RESP = charge['Restrained Electrostatic Potential Charges']
    if fname == None:
        fname = 'mol'
    else:
        fname = fname
    Chem.MolToPDBFile(pk.mol, '{}.pdb'.format(fname))
    pdbobj = PandasPdb().read_pdb('{}.pdb'.format(fname))
   # Following step is editing state, you can see the approach like a pandas method ;)
    pdbobj.df['HETATM']['b_factor'] = RESP
    pdbobj.to_pdb('{}.pdb'.format(fname))

Run the function! After running the code, I could get two files named aceticacid.pdb and tetrazole.pdb

It can load from pymol.

make_pdb_with_Resp('CC(=O)O','acetic_acid')
make_pdb_with_Resp('C1=NNN=N1','tetrazole')

Now I would like to color the molecule with edited b_factor. So I use spectrum command on pymol command line.

spectrum b, blue_red, minimum=-1., maximum=1.
Acetic acid with RESP charge

Acetic acid has negative charge on two oxygen atoms. And carbonyl carbon has the most positive charge.

Tetrazole with RESP charge

The result of tetrazole is strange for me because the nitrogen at position 2 which has hydrogen shows positive. And the nitrogen at position 1 is the most negative. I investigated two protonated state of the tetrazole, 1H-tetrazole and 2H-tetrazole.
But both results shows the nitrogen atom which has hydrogen has positive charge than other nitrogen atoms. I found that BioPandas is useful. But I could not get solution for the last problem. Hmm…

Advertisements

Draw HOMO LUMO with psikit-2 #chemoinformatics

Yesterday I wrote a post about psikit function for HOMO LUMO drawing. And @fmkz__ gave me very suggestive comment.

He performed QM calculation about tetrazole and disclosed the its distribution of negative charge.

I had interested in his suggestion and I tried it. By using psikit, I calculated energy of acetic acid and 5-methyl-1H-tetrazole and got HOMO LUMO cube files from the results.

Results is below. At first HOMO and LUMO of acetic acid.

HOMO of acetic acid
LUMO of acetic acid

These results shows that carbonyl C=O has large MO and it means carbonyl is most reactive position of the acetic acid. And also LUMO around the alpha carbon is large. It indicates that the site is reactive with electrophile I think.

Next, let’s see result of tetrazole. Tetrazole is often used in medicinal chemistry as a bioisoster of carboxylic acid.

HOMO of tetrazole
LUMO of tetrazole

These results shows 3-N position has large MO compared to other nitrogen atoms. I think it indicates that 3-N position is more reactive than other nitrogen. And I found good example to explain the result in slide share. URL is below.

https://www.slideshare.net/ShrikantPharande/tetrazole-and-triazole-67637231
Page 21 of the slide shows example of 5-sub tetrazoles alkylation.
The page shows example of reaction between 5-Me-1H-tetrazole and trityl-Cl. And most of alkylation occurred at n-3 position.

Borrowed from slide share

It is very exciting for me because I could know that QM is very useful tool for reactivity prediction. The result motivates me to enhance of psikit. ;-)

Any suggestion and advices are highly appreciated.

Draw HOMO LUMO with psikit #RDKit #Psi4 #PyMol

Now I and @fmkz___ are developing a thin wrapper library for Psi4 and RDKit named psikit.
Today I added new function for viewing molecular orbital HOMO and LUMO with pymol. Psi4 has the function which can output cube format file of MO named ‘cubeprop’. So I try to implement the function in to psikit.
By using the library, you can get HOMO/LUMO cube file easily.
Example is below.

from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from psikit import Psikit
smi = 'COc1cccnc1'
mol = Chem.MolFromSmiles(smi)
AllChem.Compute2DCoords(mol)
mol

Just call ‘getMOview’ after optimization. It took a little bit long time on my PC.

pk = Psikit()
pk.read_from_smiles(smi)
%time pk.optimize()
>Optimizer: Optimization complete!
>CPU times: user 1min 54s, sys: 3.97 s, total: 1min 58s
>Wall time: 31.5 s
>-360.58490537900167
# get MO view!
pk.getMOview()

After executing the method, psikit makes 5 files named target.mol, Psi_a_n_n-A.cube, Psi_b_n_n-A.cube, Psi_a_n+1_n+1-A.cube, Psi_b_n+1_n+1-A.cube. First file is mol file. Second and third is information of its HOMO and forth, fifth is information of LUMO.

Then load all files from pymol. Change HOMO/LUMO object show setting as dot. I could get following image.

Optimized molecule
Molecule and HOMO
Molecule and LUMO
Molecule and HOMO/LUMO

The function and pymol seem to work fine. But I want to compare data between this result and another QM’s result.
Current version of psikit is still under development and some bugs. I would like to fix and tune the code.