Hi Stamatia,
Assuming that you are interested in getting the energies of non-minimized
conformers, here's a relatively easy way to do what you want. Notice that
I'm using the MMFF94 force field, since it's probably a better choice for
energies. You could do the same thing with UFF if you'd prefer.
In [4]: m = Chem.AddHs(Chem.MolFromSmiles('OCCCN'))
In [5]: cids = AllChem.EmbedMultipleConfs(m,10)
In [6]: ps = AllChem.MMFFGetMoleculeProperties(m)
In [8]: for cid in cids:
...: ff = AllChem.MMFFGetMoleculeForceField(m,ps,confId=cid)
...: ff.Initialize()
...: print(ff.CalcEnergy())
...:
...:
12.233357032097315
4.7643146024595335
6.163107627780452
5.476964847773123
4.412833444203187
5.665311226271907
7.536117466355753
3.421494301728207
8.97909661823493
10.863891606189476
Your first approach is inefficient and just a bit wrong: the call to
MolFromMolBlock() will remove the Hs from your molecule. If you want to go
that route anyway, you should do:
b.append(Chem.MolFromMolBlock(a,removeHs=False))
Note that if you also want to minimize the conformers, you can do this:
In [9]: AllChem.MMFFOptimizeMoleculeConfs(m)
Out[9]:
[(0, 1.0966514168737902),
(0, -1.512072483231206),
(0, 0.6847373872800766),
(0, -0.6067431133406469),
(0, -1.5142123177302143),
(0, -2.08974298611677),
(0, -2.089742985761273),
(0, -2.053178209868828),
(0, -3.9515977883892077),
(0, -1.0767726284668093)]
There you get a 2-tuple for each conformer. The first element lets you know
whether or not the minimization converged, and the second is the energy of
the conformer.
I hope this helps,
-greg
On Tue, Feb 12, 2019 at 5:05 AM Stamatia Zavitsanou <
stamatia.zavitsa...@oriel.ox.ac.uk> wrote:
> Hello everyone,
>
>
> I am a bit confused. I want to make conformers and calculate their
> energies, but I am not very sure which is the correct way and at which
> point I need to embed my molecule.
>
>
>
> I am doing this :
>
>
> for line2 in infile2:
>
>
> mol2 = Chem.MolFromSmiles(line2)
>
> m2 = Chem.AddHs(mol2)
>
>
> cids = AllChem.EmbedMultipleConfs(m2, numConfs)
>
> for k in cids:
>
> a = Chem.MolToMolBlock(m2, confId=k)
>
> b.append(Chem.MolFromMolBlock(a))
>
>
> fff = AllChem.UFFGetMoleculeForceField(b[k])
>
> fff.Initialize()
>
> fff.Minimize()
>
>
> energy_value_U = fff.CalcEnergy()
>
> print("for each conformer")
>
> print(energy_value_U)
>
> energies.append(float(energy_value_U))
>
>
> Should I be doing something else? like:
>
>
> for line2 in infile2:
>
>
> mol2 = Chem.MolFromSmiles(line2)
>
> m2 = Chem.AddHs(mol2)
>
>
> cids = AllChem.EmbedMultipleConfs(m2, numConfs)
>
> for k in cids:
>
> a = Chem.MolToMolBlock(m2, confId=k)
>
> b.append(Chem.MolFromMolBlock(a))
>
>
> AllChem.EmbedMolecule(b[k])
>
> fff = AllChem.UFFGetMoleculeForceField(b[k])
>
> fff.Initialize()
>
> fff.Minimize()
>
> energy_value_U = fff.CalcEnergy()
>
> print("for each conformer")
>
> print(energy_value_U)
>
> energies.append(float(energy_value_U))
>
>
> or like this:
>
>
> for line2 in infile2:
>
> print(line2)
>
> print(i)
>
> mol2 = Chem.MolFromSmiles(line2)
>
> m2 = Chem.AddHs(mol2)
>
> AllChem.EmbedMolecule(m2)
>
> cids = AllChem.EmbedMultipleConfs(m2, numConfs)
>
> for k in cids:
>
> a = Chem.MolToMolBlock(m2, confId=k)
>
> b.append(Chem.MolFromMolBlock(a))
>
>
> AllChem.EmbedMolecule(b[k])
>
> fff = AllChem.UFFGetMoleculeForceField(b[k])
>
> fff.Initialize()
>
> fff.Minimize()
>
> energy_value_U = fff.CalcEnergy()
>
> print("for each conformer")
>
> print(energy_value_U)
>
> energies.append(float(energy_value_U))
>
>
> or maybe non of the above??
>
>
> Thanks a lot,
>
> Stamatia Zavitsanou
> _______________________________________________
> Rdkit-discuss mailing list
> Rdkit-discuss@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
>
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
