failing to read PDB files generated by VMD

[lmiq]

When trying to read this file (and other similar ones), I get:

julia> read("vmd.pdb", BioStructures.PDB)
ERROR: Two copies of the same atom have the same alternative location ID. Existing atom:
Atom N with serial 1, coordinates [-41.156, -40.019, 21.411]
New atom record to add:
AtomRecord N with serial 1501, coordinates [-42.523, -41.722, -19.881]

This is one example VMD pdb file.

vmd.pdb.zip

[jgreener64]

The offending PDB lines are

ATOM      1  N   POPE    1     -41.156 -40.019  21.411  1.00  0.00      L11  N

and

ATOM   1501  N   POPE    1     -42.523 -41.722 -19.881  1.00  0.00      L21  N

which both try to put an atom with name N onto residue 1 of the chain with empty chain ID.

We treat this defensively since there are no or very few cases in the PDB that have such duplicate atoms. I reported a handful of cases to the PDB a while ago and they updated the records, so I believe it is considered a format violation.

In this case you could use different chain IDs (the empty column 22) or keep incrementing the residue number.

[timholy]

Since this file was written by VMD, I assume this should be reported to them? @lmiq, can you do that?

[lmiq]

Well, there is nothing wrong with that file for the use it was designed for, so I do not think there's anything to report there, upstream. The parsing of the residues in that file is simply done by incrementing the residue counter when the residue index changes (for more or less). Having limitations in the number of residues per chain, or number of chains, etc, is something that cannot be important in MD simulations PDB files.

This is the choice to be made here: have or not the possibility of parsing non-standard files to some degree.

I perfectly understand if BioStructures keeps strictly adhering to the standard, but that would not allow us using it in the most broad context of MD simulations.

Note that this limitation, specifically, is associated to trying to read the data into the hierarchical structure, so in some sense this is related to that initial choice of representation.

Also, note the L11 and L21 segment identifiers. These are used as the "top" level of hierarchy in the context of VMD and some simulations packages, but are barely used for reporting experimental structures.

[timholy]

If it's true that it violates https://www.wwpdb.org/documentation/file-format (I haven't checked), I'd say that's clearly a problem. It's not OK that there's a workaround. It introduces ambiguities and issues, precisely as being discussed here.

But if the format is ambiguous, then that's another issue entirely.

[jgreener64]

It's certainly challenging to balance the considerations.

Note that this limitation, specifically, is associated to trying to read the data into the hierarchical structure, so in some sense this is related to that initial choice of representation.

This is the key, BioStructures guarantees that every atom has "meaning" about what it is (i.e. on which residue and chain), which is useful in many contexts (such as file interchange) but in turn necessitates some complexities in representation and strictness in parsing.

The only way the original case could work currently is if the existing atom is overwritten (or the new atom ignored) with a warning, which is unlikely to be desired behaviour even if it runs without error.

I guess philosophically the aim of BioStructures is to not so much to read structural files, but to represent unambiguously the molecules within them. This affects other design decisions, such as why countatoms will treat disordered versions of the same atom as one (unless expand_disordered is used).

But if the format is ambiguous, then that's another issue entirely.

Sadly, this seems to be the case. There is a lot of documentation on the column format, but as far as I know not so much on the row format (e.g. what duplicates mean). In general I prefer it when other tool authors use MODEL blocks to distinguish conformations in a MD trajectory though. That works well in BioStructures.

[lmiq]

From the description, that file does not adhere to the standard because:

- Non-blank alphanumerical character is used for chain identifier.

and there the chain identifier is blank. But the standard does not say anything about two residues in the same chain having the same number (which is actually the issue here). For actual bio-structures it is just assumed that that does not make sense. But when part of the system is thousands of water molecules, that is just unnecessarily limiting.

[timholy]

Given that the format is ambiguous, then perhaps we do need to deal with ambiguity. I don't think this should be open-ended: it's kind of a disaster to take the strategy "can I find something that seems to make sense in column 5? No? OK, let's see if things seem more sensible if I substitute the value in column 11? No? OK, maybe I can infer a good value from column 13?" But if instead one passes flavor=:vmd to the reader, and that follows strict rules, things are more sensible.

[timholy]

when part of the system is thousands of water molecules

is this part of the issue? That (legacy) PDB only lets you use a single character for encoding the chain? In which case shouldn't VMD be writing mmCIF instead?

[lmiq]

As a minimal example:

ATOM   4247  H   HOH     1     -28.310   3.370 -27.343  1.00  0.00      HATI                          
ATOM   4248  H   HOH     1     -27.809   4.855 -27.803  0.00  0.00      HATI                          
ATOM   4249  O   HOH     1     -27.550   4.020 -27.318  0.00  0.00      HATI                          
ATOM   4250  H   HOH     2     -19.010   5.493  44.300  0.00  0.00      HATI                          
ATOM   4251  H   HOH     2     -18.958   5.838  45.896  1.00  0.00      HATI                          
ATOM   4252  O   HOH     2     -19.119   6.215  44.984  0.00  0.00      HATI                          
ATOM   4253  H   HOH     1      12.951   6.641  28.351  1.00  0.00      HATI                          
ATOM   4254  H   HOH     1      13.590   7.926  29.130  0.00  0.00      HATI                          
ATOM   4255  O   HOH     1      13.304   6.973  29.226  0.00  0.00      HATI 

VMD, Pymol, Packmol, and other common simulation packages interpret this without problems. There are three water molecules there.

The issue is that this does not fit into an hierarchical structure, unless the reader creates chains arbitrarily. That would come with a bunch of problems as well.

is this part of the issue? That (legacy) PDB only lets you use a single character for encoding the chain? In which case shouldn't VMD be writing mmCIF instead?

Well, in some sense yes, the PDB format is problematic. But that ship has sailed decades ago in the MD field. Not that there aren't other formats that are more appropriate, particularly for the limitations of the coordinates fields, but the legacy PDB format is still widely used because of its readability. This is something we (MD simulation package authors) just have to live with.

ps: MIToS is reading this if one adds occupancy and b-factor fields (added now), because its hierarchy is at the residue level only, and it just not cares about the "meaning" of two residues with identical numbers in the same chain.

[timholy]

So maybe we need a flavor=:md mode or something? @jgreener64 would that be acceptable?

[jgreener64]

Yes, but the question still remains about how to put the two atoms into the same hierarchical structure.

We could assign an unused residue number or chain ID, but that feels clunky and runs into problems if the assigned value is found later during reading. We could have another data type that stores multiple coordinates for the same atom, but that raises concerns about how to write the file back out in the same order. We could assign duplicates into a new Model, which would mimic MODEL/ENDMDL tags but doesn't fit the case of many waters in a MD file.

[lmiq]

Maybe the flavor should be the type of hierarchy desired: :segment, :chain, :residue (MIToS now), :atom (PDBTools) or whatever. I don't know, however, if that would make most of the other functionality of BioStructures fail to interoperate now, and thus the package very complex.

[timholy]

I guess my thought is that we decide on concrete flavor-specific behavior and document it. Passing a flavor might require that the entire file be scanned first to identify "occupied" names/indices/whatever, which would be a bit of a performance hit, but perhaps the price of abusing a file format.

[lmiq]

I think it is a little more complicated than that. If we keep the hierarchy we are forced to attribute some of the fields arbitrarily. That will cause confusion.

[timholy]

Isn't the hierarchy about how you assign bonds? Meaning, you have implicit hierarchy any time you want to link entries in rows together.