Amber mapping ignores hydrogen in BB

[fgrunewald]

The Bug

It seems the amber mapping files do not produce the expected mapped coordinates. The difference is typically very small (i.e. sub Angstrom), which is probably why it wasn't noticed. However, this might cause problems later down the line especially when using the martinized PDB file as reference for the EN other other things.

Files to reproduce the bug

Simple ALA input PDB.

TITLE     PSI
REMARK    THIS IS A SIMULATION BOX
CRYST1  101.000  101.000  101.000  90.00  90.00  90.00 P 1           1
MODEL        1
ATOM      1  N    ALA  652      15.180   2.910   2.110  1.00  0.00           N
ATOM      2  H    ALA  652      15.030   3.680   2.730  1.00  0.00           H
ATOM      3  CA   ALA  652      14.680   1.590   2.480  1.00  0.00           C
ATOM      4  HA   ALA  652      15.450   0.960   2.340  1.00  0.00           H
ATOM      5  CB   ALA  652      14.270   1.570   3.950  1.00  0.00           C
ATOM      6  HB1  ALA  652      13.930   0.670   4.190  1.00  0.00           H
ATOM      7  HB2  ALA  652      15.050   1.800   4.520  1.00  0.00           H
ATOM      8  HB3  ALA  652      13.540   2.250   4.090  1.00  0.00           H
ATOM      9  C    ALA  652      13.520   1.160   1.590  1.00  0.00           C
ATOM     10  O    ALA  652      13.690   0.380   0.660  1.00  0.00           O
TER
ENDMDL

This should map to:

ATOM      1  SC1 ALA X   1      14.198   1.572   4.188  1.00  0.00              
ATOM      2  BB  ASP X   2       11.397   1.038   1.498  1.00  0.00              

But it actually maps to:

ATOM      1 BB   ALA     1      14.259   1.478   1.647  1.00  0.00              
ATOM      2 SC1  ALA     1      14.270   1.570   3.950  1.00  0.00              

The relevant command is:

martinize2 -f ALA.pdb -ff martini3001 -x cg_ala.pdb -from amber

Solution

Remove the '!' from the mapping files. However, this seems to be in contradiction to what their actual purpose was (i.e. indicate optional hydrogen atoms to map). As I understand it optional means use them if present. @pckroon did I understand the mapping file syntax incorrectly.

[pckroon]

No. The "!" means give these 0 weight. There is no such thing as optional in the mappings. And it doesn't have to be, since the graph is always complete.

[fgrunewald]

I actually moved this to the discussions board, because I'm not 100% if this is a bug or design feature. It seems a bug to me because the side-chain bead positions would be incorrect. @paulocts can you comment, if the hydrogen atoms should be omitted.

[paulocts]

I would be fine with positions slightly different if this feature would be corrected/reduced in the minimization step, which is the case for all bonded terms, with three possible exceptions: elastic networks, go models and side chain corrections. As this choice of including H or not affects the equilibrium values and, in the case of elastic/Go even the presence of the bond/LJ potential, then I would say this is not good. 

This behaviour would mean that the model would be different because of the differences in the positions of the beads coming from the presence/absence of hydrogens. In this aspect, I guess the safest and easiest options would be two excludent possibilities:

1. the structures always need to have all the hydrogens there for martinize2

or

2. we never consider the hydrogen atoms for the mapping in martinize2

The option 1 may help to be sure that the models can always be minimized as the the positions of the beads will reflect more the bond lengths, angles, etc of the model as nowadays the Martini 3 models include the hydrogens to better reproduce SASA and molecules volume. 

The ption 2 is the one that makes it easier for the users as they would not need to add hydrogens in case the pdb is already good. I think is also more consistent with how the models were generated in the olde martinize for Martini 2. You may need a bit more minimization steps to relax the structure, but I think it works fine.

The important thing is to be able to reproduce the models/data which now seems compromised. Essentially, if the users do not informed in their method sections of the paper if they include H or not, other people may not be able to reproduce the model.

PS: There is a option 3 which would be to make martinize2 to add H atoms in case they are missing. But this would demand too much extra coding from your side.

[pckroon]

I'm going to veto option 1. It will make M2 too hard to use. I disagree with the statement that it blocks reproducibility. User should always report their initial atomistic input structure anyway, and that won't magically gain or lose hydrogens. (and there's also the -ignh flag)

I do agree that it's an issue though. I think it would be a good idea (in general) to check whether any NaN positions are used to generate bead positions, and warn if those positions get used to determine equilibrium values. Note that this does not interfere with option 3, since it may also apply to e.g. missing side chains.

[fgrunewald]

So the choice is unproblematic in relation to unbiased proteins. That's why I moved this to the discussion. However, it becomes a problem if the EN is added, because the mapped EN is not consistent with the equilibrium values of the backbone. That introduces some strain due to competition between EN and BB bonded terms. For larger proteins this becomes noticeable. I bet this leads to some of the instabilities we see. Energy minimisation cannot fix this.

In my opinion, for Martini3 the -ignh flag should raise a warning along the lines of You are missing hydrogen atoms this will effect your structural biases like EN and lead to a starting structure that is not consistent with Martini3. That would be easy to implement.

I think option 1 light might be the best way to go. We simply require or give the backbone hydrogens a weight of 1. HN is already required, so it would be easy enough to add HA. That would solve most issues in relation to EN and structural biases. The only way to overrule this choice would then be to add -ignh but that would raise a warning.

The only remaining issue is that I've seen people using martinize2 to forward map structures in oder to convert trajectories to Martini resolution. Luis did that in the PTF I think, but also other people are doing this. Of course that will lead to structures which are not equivalent to the actual trajectories of a simulation. I'm not sure how to handle this best though.

@pckroon can we raise a warning somewhere IF hydrogens are missing in combination with using Martini3?

[pckroon]

The problem is broader: consider a protein in which backbone atoms are missing, and you need to construct the EN with martini2.
Issuing a warning for the combination martini3 + -ignh + EN is a good idea. But I still think that the general solution is to issue a warning for all cases where coordinates determine parameters, and some bead coordinates were constructed from NaNs.
Just adding atoms to mappings is not enough, since we reconstruct all graph nodes anyway.

I'm not sure what you mean with the trajectory mapping. That's not really a supported use case anyway though.

[fgrunewald]

You mean constructed from NaNs at the AA level? Because the real problem is when atoms are missing but there are still bead coordinates generated at the CG level and the structural bias is based on those.

By trajectory mapping I mean, that people literally write a trajectory as PDB and use martinize2 to coarse-grain them, because there currently (with the exception of the less known fast-forward) no program for doing so.

[pckroon]

You mean constructed from NaNs at the AA level? Because the real problem is when atoms are missing but there are still bead coordinates generated at the CG level and the structural bias is based on those.

This is the same thing? Atoms that are missing at the AA level get constructed with NaN coordinates, and those thus contribute to the bead coordinates. Fortunately not all bead coordinates affect the topology.

By trajectory mapping I mean, that people literally write a trajectory as PDB and use martinize2 to coarse-grain them, because there currently (with the exception of the less known fast-forward) no program for doing so.

Ouch. Do you know anyone who wants an easy paper? Make mapping matrix once, have MDAnalysis do a dotproduct.