I think the following will work:
for i in range(len(results.ligands)):
ligand = results.ligands[i]
I'm curious as to why the docking results are corrupted, though. Did you concatenate many files to make the docking result, or did GOLD create corrupted files?
can you let me know what the error message is when you try to access the occupancy attribute, and give me details of your Manjaro operating system? I'd like to investigate this if I can find an appropriate vagrant virtualbox for it.
yes you will be able to use the API for this.
The docking module of the CSD API, https://downloads.ccdc.cam.ac.uk/documentation/API/descriptive_docs/docking.html
shows how to read ligands in their docked poses, how to get access to the scoring terms of the fitness function you have used, and how to inspect the hydrogen bonds of the docked poses.
There are example programs using docking in https://downloads.ccdc.cam.ac.uk/documentation/API/cookbook_examples/docking_examples.html which may be helpful.
Can you show me an example script where atom.occupancy gives you an error? The occupancy attribute may be None if the atom has no coordinates, that is for an implicit hydrogen.
I can use it perfectly normally in, for example:
from ccdc import io
csd = io.EntryReader('csd')
mol = csd.molecule('ABABUB')
print('\n'.join('%s: %.4f' % (a, a.occupancy) for a in mol.atoms))
The value 'Mixed' is returned when the algorithm used to determine chirality is unable to come up with an answer. This can occur when the chiral atom has more than four bonds, for example, atom U1 of AAPUNI.
the entry for csd number 1827879 will be appearing in the February update to the CSD which should be available to you very soon.
what problems are you having with the script?
the return value of the packing() function on a Crystal is an instance of a Molecule, not a Crystal, and the 'intermolecular' attribute of the HBondCriterion does not apply: all the HBonds are considered to be intramolecular, even if they are in disconnected components of the molecule. Additionally, the HBonds will not be found since the path_length_range attribute if the HBondCriterion is left at its default setting of (4, 999), requiring the HBond to be within a single connected component of the molecule.
You can discover the HBonds by setting the path_length_range to (0, 999), allowing HBonds to be between disconnected components of the molecule, and you can determine which HBonds are between disconnected components by testing the shortest_path() method of the atoms of the HBond, that is to say
criterion.path_length_range = (0, 999)
criterion.intermolecular = 'ANY'
packed = cry.packing()
hbonds = packed.hbonds(hbond_criterion=criterion)
inter_component_hbonds = [hb for hb in hbonds if packed.shortest_path(hb.atoms, hb.atoms[-1]) == 0]
I hope this is helpful, If anything is unclear please get back to us.
the definition of what constitutes an hbond is contained in the class molecule.MoleculeHBondCriterion. Here there are options for the geometry, the presence or absence of the hydrogen, and crucially for you, a means to specify atom types which may form a hydrogen bond.
These are documented in the API in:
I hope this is helpful.
you can use the python API to do this. For example:
from ccdc import io
with io.EntryReader('/path/to/file.sdf') as reader, io.EntryWriter('/path/to/file.pdb') as writer:
for entry in reader:
I hope this is helpful. Please come back to me if anything needs clarification.
I don't have an existing script to do this, but here is a sketch of how you might go about it.
from ccdc import io, descriptors
directory = '/path/to/directory/of/cifs'
cif_files = glob.glob(os.path.join(directory, '*.cif'))
for cif_file in cif_files:
with io.MoleculeReader(cif_file) as rdr:
cif = rdr
for a1 in cif.atoms:
for a2 in cif.atoms:
print a1.label, a2.label, descriptors.MolecularDescriptors.atom_distance(a1, a2)
It is worth noting that if you use an EntryReader to read the CIF file, all of the CIF attributes are available, for example,
cif_entry = io.EntryReader(cif_file)
a1 = cif_entry.attributes['_geom_bond_atom_site_label_1']
a2 = cif_entry.attributes['_geom_bond_atom_site_label_2']
length = cif_entry.attributes['_geom_bond_distance']
for a1, a2, length in zip(a1, a2, length):
print a1, a2, length
Hope this is helpful. Please come back with any questions you have about this, or indeed any other API matter.