I am relaxing a dimer structure which contains one modified residue per subunit. Relaxing a monomer with one modified residue is OK, as is relaxing a dimer where the modified residue has been removed from the input PDB file. However, as soon as there are two copies of the modified residue in the input PDB or during relax run (a dimer or a monomer with C2 symmetry enabled) then the calculation run into this dead end:
ERROR: dis==0 in pairtermderiv!
ERROR:: Exit from: src/core/scoring/methods/PairEnergy.cc line: 468
Reading the offending code line it appears that an interatomic vector is sensed as having a length of zero hence the error. To the best of my knowledge, the input dimer structure I am using, even with H added (using reduce -BUILD), has no interatomic distance that short nor any other serious clash. I also checked there were no duplicated residues.
Since the problem goes away as soon as I remove the modified residue, I am guessing that one is the problem... even though it gets relaxed when there is a single copy of it during the run.
I am stuck at this point, so any help will be much appreciated!
PS BTW the above error kicks in right after the following output lines:
protocols.relax.FastRelax: CMD: repeat 8328.73 0 0 0.44
core.pack.interaction_graph.interaction_graph_factory: Instantiating DensePDInteractionGraph
core.pack.pack_rotamers: built 30895 rotamers at 743 positions.
core.pack.pack_rotamers: IG: 103244808 bytes
This suggests this problem might be connected to packing rotamers (I have not defined rotamers for my modified residue) but again, the correct termination of the relax job using one copy of the modified residue would speak against that...
"ERROR: dis==0 in pairtermderiv!"
I am guessing you are running into an issue where the fa_pair term cannot deal with the pair of your two modified residues. With only one residue present, all pairs are normal-normal or normal-modified; with the modified residue present then modified-modified can occur, and for some reason it comes up bad.
The pair term acts on the "ACTCOORDS" of a residue. Is there a record field for that in your residue type? I'm not sure exactly what the field is called - try comparing it to the unmodified residue type.
When I diff the LYS.params file against my modified residue params file I see that the ACT_COORD_ATOMS field is actually missing in the latter:
LYS.params => ACT_COORD_ATOMS NZ END
LLP.params => (missing)
Nearby lines are all equal except for this one, whose relevance escapes me for now:
LYS.params => NBR_RADIUS 5.0084
LLP.params => NBR_RADIUS 16.572544
I am going to try to add the ACT_COORD_ATOMS line into my modified residue's params file and re-run these jobs, posting the result back to the forum.
Thanks a lot for the pointers!
So long as NZ exists, feel free to copy the ACT_COORD_ATOMS line. Otherwise use whatever atom you feel like (maybe CB). The pair term should be returning zero for a modified residue anyway (no pair statistics to calculate from...).
NBR_RADIUS is used in calculating interaction graphs.
Most two body energy terms in Rosetta have a certain distance cutoff - a distance outside of which the value drops to zero (this is what we mean when we say "short range energy term"). We take advantage of this in certain cases so that we don't have to compute all the residue-residue interactions when rescoring - just the ones that are close enough. A trick to this is that you need to know how far apart each atom is. That's where the NBR_RADIUS comes in. It's the furthest distance from the neighbor atom that a heavy atom will be in any of the available rotamers. Because the neighbor atom does not move in packing, you know that two residues with neighbor atoms greater than (NBR_RADIUS_A + max_interaction_dist + NBR_RADIUS_B) apart from each other will always have zero short-range two body energies, so you don't have to waste time recomputing them.
That's why the NBR_RADIUS is bigger on your modified lysine - there's more stuff hanging off of it, so heavy atoms can reach further away from the neighbor atom.
OK, so I did exactly that (copying the ACT_COORD_ATOMS field from lysine over to my modified residue's params file) and started the relax job without changing the scripts nor the input files, and this time the job has being running for about 3 min without stopping, which is a good sign that it has finally taken off (previously it was grinding to a halt within a few seconds).
Thanks so much!
This is simply to confirm that the above solution works correctly and produces the intended output at the end of the job.