chryswoods
commented
1 year ago Hi - thanks for your kind words :-)
The amber.zmatrices file is a z-matrix file. This file doesn't specify a forcefield. Instead, it tells waterswap how to perform the Monte Carlo moves that are used to sample the system.
Specifically, it tells waterswap how to perform the Monte Carlo moves on the protein molecules. (it doesn't need to know how to move the waters, as these are only translated and rotated rigidly, and it uses an algorithm to automatically generate the Monte Carlo moves for small molecules like ligands or cofactors - these algorithms don't work for proteins).
Every protein residue is divided into a backbone (containing bbatoms) and sidechain (containing scatoms). The sidechains are the same regardless of the position of the residue in the protein (e.g. C-terminus, N-terminus or middle).
The backbone atoms - which are the C, N, O and connected hydrogens of the backbone - are specified as "chain" templates, e.g. chain aamiddle describes how general amino acid backbone atoms would move for residues in the middle of the protein, while chain aanterm does the same for backbone atoms in the residue at the N-terminus.
The residue templates, e.g. residue ALA say how to move the sidechain atoms of the residues. For example, residue ALA gives the names of the chain templates used for its backbone atoms at different chain positions, then it gives the z-matrix that defines how each of the sidechain atoms are built and moved from those backbone atoms.
The error message shows that an atom called H exists in one of your protein residues which isn't expected based on the atoms listed for that residues residue and chain (backbone) template.
Typically I have seen these errors if the protein residues don't use standard (PDB-style) atom names, or if the groups on the C- or N- terminals of the protein aren't specified in the amber.zmatrices file.
The question to answer are;
-
What are the terminal residues of your protein? Only NME, ACE or standard amino acids are current specified in the
amber.zmatricesfile -
How are the hydrogen atoms named in your protein? You need to match the names in the templates, e.g.
residue ALA
rigidbody rotate 2.000 translate 0.050
bbatom N
bbatom CA
bbatom O
bbatom C
backbone first aanterm middle aamiddle last aacterm single aasingle
zmatrix CB CA N C
zmatrix HB1 CB CA N
zmatrix HB2 CB CA HB1
zmatrix HB3 CB CA HB1
angle CB CA N flex 0.50
angle HB1 CB CA flex 0.50
angle HB2 CB CA flex 0.50
angle HB3 CB CA flex 0.50
dihedral HB1 CB CA N flex 10.00shows that the sidechain hydrogens in ALA need to be called HB1, HB2 and HB3.
For terminal residues, the naming also needs to be precise, e.g. here's the templates for NME
chain nme
zmatrix HH31 CH3 N H
angle HH31 CH3 N flex 0.50
dihedral HH31 CH3 N H flex 15.0
residue NME
rigidbody rotate 2.000 translate 0.050
bbatom H
bbatom N
bbatom CH3
bbatom HH31
backbone first nme middle nme last nme single nme
zmatrix HH32 CH3 N HH31
zmatrix HH33 CH3 N HH31
angle HH32 CH3 N flex 0.50
angle HH33 CH3 N flex 0.50showing that the hydrogens need to be called H, HH31, HH32 and HH33.
I'm sorry that waterswap isn't more helpful in its error message. Unfortunately there won't be any more useful details, even if get_last_error_details() worked for the old API. This is because those details give a backtrace and thread information, which historically confused more people than they helped (hence why they are now hidden). The backtrace won't give any more info about which residue has the wrong atom.
The way to debug this is to look for the H atoms in your protein and make sure that they appear in the chain and residue template for the associated residue. Start with the terminal residues first, as these are normally the source of the problem.
I will say though that waterswap is quite an old program (2013). It is based on Monte Carlo and so can't be GPU-acclerated. It will likely be a lot slower than a GPU-accelerated MD code and free energy method, and the free energies will not be as high a quality as a more modern metadynamics or other biased dynamics method (e.g. you may have more success with alchemical transfer methods, e.g. https://github.com/Gallicchio-Lab/AToM-OpenMM).
ljmartin
Hi all - thanks for making this great software available.
I'd like to try out waterswap, and have installed it in a container using the instructions here. It seems to run happily on example inputs given in a tutorial (here), but on my own inputs, generated via OpenMM and OpenMMForceFields, it can't generate a Z matrix for something.
I guess this is because my residues don't follow the naming scheme for the file amber.zmatrices. However, I can't figure out which residue is incorrect. The error is:
I can guess is that there's an atom or a residue called
Hthat shouldn't be. There aren't any residues calledH(I checked with parmed), and atomHis, I think, a backbone amide hydrogen atom so shouldn't be a problem ordinarily.I can reach this stage in an IPython shell, but I can't
call sire.error.get_last_error_details() for more infobecause Waterswap runs within the old API and this:sire.error.must be in the new API, since it complains about trying to run modules from two different APIs.Zmatrixmaker is in C++ so I can't 'debug' further. Is there any other way I can track down what the residue is? Also, what forcefield is amber.zmatrices? TY!