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 has negative charge on two oxygen atoms. And carbonyl carbon has the most positive 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…
