I have a list of SMILES of small molecules and I want to be able to simulate these molecules with an Amber force field in Python. Currently, I use RDkit to convert the smile into a PDB file:
from rdkit import Chem
from rdkit.Chem import AllChem
molecule = Chem.MolFromSMiles('smile')
molecule = Chem.AddHs(molecule)
AllChem.EmbedMolecule(molecule)
Chem.MolToPDBFile(molecule, 'molecule.pdb')
Then I use the ParmEd package:
import parmed as pmd
pdb = pmd.load_file('molecule.pdb')
This pdb object for 'CC(C)C' for example, looks like:
<Structure 14 atoms; 1 residues; 13 bonds; NOT parameterized>
Then I want to create a prmtop file, which is needed to initialise an Amber force:
pdb.save('molecule.prmtop', overwrite=True)
This throws the error:
RuntimeError: atom types are not assigned
I'm not sure why this error gets thrown, since the atom types do seem to be inside the PDB file. I'm new to trying to do these kinds of things, so I'm not completely sure if more information is needed.
If this method in general is not preferable, what is the best way to achieve this (preferably on python)?
Any help would be greatly appreciated.
