Erroneous interplay between gmx rms command and atommass.dat: Can not find mass in database for atom MG in residue - Redmine #3368
Archive from user: Vedat Durmaz
Bug is related to the closed Gromacs Bug #20 (issued by Erik Lindahl on
10/18/2005 04:44 PM):
“g_rms uses incorrect masses when reference is a PDB/GRO file”
The following description was copied from a mail to the Gromacs users list:
I’m pretty sure it’s not a feature but a bug which I’ve faced in the GMX versions 2018.7, 2019.5 and 2020.
When I try to calculate RMSD values for a protein system including a catalytic magnesium ion “Mg” using the command
$ gmx rms -s mol.pdb -f mol.xtc -f2 mol.xtc -o mol-rmsd.xvg -debug
i get this error:
<code class="text">
can not find mass in database for atom mg in residue 1370 mg
-------------------------------------------------------
program: gmx rms, version 2019.5
source file: src/gromacs/fileio/confio.cpp (line 517)
fatal error:
masses were requested, but for some atom(s) masses could not be found in the
database. use a tpr file as input, if possible, or add these atoms to the mass
database.
let’s accept for the moment that using a .tpr file with -s (rather than
a .pdb file) is no option for me. consequently, gmx retrieves atom
masses from atommass.dat which actually contains the mg ion:
<code class="text">
; note: longest names match
; general atoms
; '???' or '*' matches any residue name
??? h 1.00790
??? he 4.00260
??? li 6.94100
??? be 9.01220
??? b 10.81100
??? c 12.01070
??? n 14.00670
??? o 15.99940
??? f 18.99840
??? ne 20.17970
??? na 22.98970
??? mg 24.30500 <<< in this line
??? al 26.98150
??? si 28.08550
??? p 30.97380
having a look at the log file of “gmx rms -debug …”, i see some strange output. here’s a snippet:
<code class="text">
searching residue: ??? atom: h
not successful
searching residue: ??? atom: he
match: ??? h
searching residue: ??? atom: li
not successful
searching residue: ??? atom: be
not successful
searching residue: ??? atom: b
not successful
searching residue: ??? atom: c
not successful
searching residue: ??? atom: n
not successful
searching residue: ??? atom: o
not successful
searching residue: ??? atom: f
not successful
searching residue: ??? atom: ne
match: ??? n
searching residue: ??? atom: na
match: ??? n
searching residue: ??? atom: mg
not successful
searching residue: ??? atom: al
not successful
searching residue: ??? atom: si
not successful
searching residue: ??? atom: p
not successful
searching residue: ??? atom: s
not successful
searching residue: ??? atom: cl
match: ??? c
...
searching residue: ??? atom: cu
match: ??? c
...
i don’t know what exactly the rms command is doing behind the scenes and
also don’t want to agonize over the cpp code. but to me it seems as if
the element assignment is based only on the first letter of each element
name. if i use copper (cu) instead of magnesium, the program runs fine.
now let’s compare the cu lines with the mg lines in the debugging output
above.
<code class="text">
searching residue: ??? atom: cu
match: ??? c
<code class="text">
searching residue: ??? atom: mg
not successful
obviously, cu is found because the first letter of its element
abbreviation, c (because cu\[index 0\] is c) does exist as an element,
but there is no element m, the first letter of mg. a little test
(respectively an improper workaround). if i add a line for the element
“m” to atommass.dat like this
<code class="text
??? Na 22.98970
??? Mg 24.30500
??? M 24.30500 <<< new line
??? Al 26.98150
then the rms command executed on the protein with Mg runs without any error. But this observation also implies that masses, the rms command retrieves from atommass.dat, are wrong because, e.g., to each of Cl, Cr, Cu and Co, the mass of C is assigned.
I see 2 options for GMX developers. Either you check the rms<->atommass.dat interplay, or you disable the possibility to use PDB files with the -s option. However, I would strongly discourage from the latter decision since there are cases where you have an xtc trajectory possibly generated with another tool along with a pdb file, but don’t want to spend too much time in the generation of a proper tpr file.
(from redmine: issue id 3368, created on 2020-02-05 by gmxdefault)