Hi @JoaoRodrigues , new formats are a great addition. Can you just create a PR and then we can comment directly on the code, run tests, etc?

Sounds good. Will do. Thanks!

Pure python implementation: https://github.com/Electrostatics/mmcif_pdbx

Also, chemfiles reads mmCIF and can be used inside MDAnalysis, see Reading trajectories with chemfiles.

Note on trying to use chemfiles:

With the 2.0.0b,

u = mda.Universe("1ake.cif.gz", topology_format="chemfiles")

fails with

'CHEMFILES' isn't a valid topology format, nor a coordinate format

I can't just read it with

u = mda.Universe("1ake.cif.gz", format="chemfiles")

either as it gives the same ValueError.

Just trying something stupid with the topology as a PDB also fails

u = mda.Universe("1ake.pdb.gz", "1ake.cif.gz", format="chemfiles")

with

ChemfilesError: can not read file '1ake.cif.gz' at step 0, it does not contain any step.

Conclusion: I couldn't get it to work with chemfiles. Maybe @Luthaf has some ideas but a key problem seems to be that our (EDIT) chemfiles converter do not work for topologies —— ??

I think that topology reading is not registered with MDA, since the chemfiles adapter is implemented as a coordinate reader/writer. Although the conversion from chemfiles to MDA topologies is already implemented, so it should be mostly a question of adding a new ChemfilesParser class. I'll have a look at this next week!

ChemfilesError: can not read file '1ake.cif.gz' at step 0, it does not contain any step.

This is a bit strange and probably a bug, is this 1ake from the wwwPDB?

This is a bit strange and probably a bug, is this 1ake from the wwwPDB?

Ok, I understand this part. The core of the issue is that to relatively different formats want to use the same extension: mmCIF and crystallography CIF. While they both use the same STAR format, they specify data in different ways. Unfortunately, chemfiles is associating the .cif extension with crystallography CIF files, and uses .mmcif for mmCIF.

There is a simple workaround though, since you can specify the format to use manually, with something like

u = mda.Universe("1ake.pdb.gz", "1ake.cif.gz", format="chemfiles", chemfiles_format="mmCIF / GZ")

Unfortunately, this still fails in the case of 1AKE (but should work for other files). I'll fix the 1AKE issue, it should be working in the next patch release.

I would also like to introduce a better format guessing functionality, to decide between mmCIF and crystallography CIF on the fly instead of having the use specify it manually.

The mmCIF/PDBx format would also be needed for alphafold #3377 .

does anyone have a solution for writing an MDAnalysis universe as PDBx/mmCIF?

you could use pdb_tocif tool from pdb-tools applied to all PDB files from a trajectory and then pdb_mkensemble to merge them.

I just tried pdb_tocif, but it screws up the segids. It should have segid of B0, or B1, ... It only keeps the chain (B). See the question marks below:

ATOM   374778 H  HA     ?     B   ? 116   ?    352.061     37.526    138.818      1.000      0.000 ? 116       B HA   1
ATOM   374779 C  CB     ?     B   ? 116   ?    353.637     37.061    140.192      1.000      0.000 ? 116       B CB   1
ATOM   374780 H  HB2    ?     B   ? 116   ?    353.967     36.114    140.621      1.000      0.000 ? 116       B HB2  1
ATOM   374781 H  HB3    ?     B   ? 116   ?    354.280     37.265    139.336      1.000      0.000 ? 116       B HB3  1
ATOM   374782 C  CG     ?     B   ? 116   ?    353.874     38.157    141.231      1.000      0.000 ? 116       B CG   1

Sorry to cross-post on a different project, but could you share an input file?pdb_tocif is not designed to handle multi-character chain identifiers. Feel free to open an issue in the pdb-tools repo.

I wonder if it'd be ok to use gemmi library (link) as a dependency to read mmcif and other cif-like formats consistently?

IMO it's one of the best supported crystallography-related libraries, is maintained by ccp4 and globalphasing, and allows very detailed cif parsing and writing.

For example, reading all atom ID and coordinates could be done as simple as:

# read only ATOM groups of chain 

model = gemmi.read_structure('6rz6.cif')[0]
residues = [res for res in model['A'] if res.het_flag == 'A'] # reading non-heteroatom atoms for simplicity
arr = np.array([[at.serial, *at.pos.tolist()] for res in residues for at in res])

ids, xyzs = arr[:, 0].astype(int), arr[:, 1:]

the whole discussion with devs here

@marinegor I think @richardjgowers has some ideas in this area

@richardjgowers could you share them (here, if it's appropriate place, or somewhere on discord)?

@marinegor I'd done this at a hackathon: #4303

The problem with this approach is that it doesn't do the "table join" on conect records that mmcif relies on, so you won't get some data (bonds).

I've since done this: https://github.com/OpenFreeEnergy/pdbinf which does do the "table join" to get bonds. It's into rdkit format, but conversion is trivial...