Suspiciably high delta H values when alanine scanning

[Paul Schrank]

Hey guys, 

I would have another question (hopefully my last). When I am running my alanine scan at certain positions it gives me very large positive free energy of binding for the mutant and delta H. In the FINAL_RESULTS.dat it additionally says: 

A/2 - ARGxGLY MUTANT:
POISSON BOLTZMANN:

WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS AND
THE VALIDITY OF THESE RESULTS ARE HIGHLY QUESTIONABLE!

Some absolute differences in the internal potential terms are greater than 0.005.
This should not happen when using Single Trajectory Protocol!

You can generate a detailed *.csv file with all the terms and differences as follows:
gmx_MMPBSA --rewrite-output -eo energy_terms_differences.csv

(the same for GB calculation).

Do you have any idea, what this could be. I found that it might have something to do with inconsistencies in the MM energy terms of the protein, when generating the mutant, is this correct. As always here is a sample job for one of the positions I carried out: 

https://gigamove.rwth-aachen.de/de/download/c64976ff71079a16e3674bb26a6022ac

Thank you and as always best regards!

Paul 

[Mario Sergio Valdes]

I've been checking this system since yesterday... It has me baffled... Everything seems to be fine because the topologies are consistent, there are no weird bonds, no unbonded atoms, and no other frequent problems. The most curious thing is that the PDB of the receptor is used to generate the topology of the receptor itself and the complex. So any problem should be reflected in both. However, this does not happen. The energy of the complex is consistent with the mutation, while that of the receptor varies greatly. I'll try to see what's going on and tell you.

[Paul Schrank]

Dear Mario,

Thank you for your support!

I tried multiple times, but it seems to no avail... But I saw something in the leap.log:

> LIG1 = loadmol2 CBZ_Asn_Gly.mol2
Loading Mol2 file: ./CBZ_Asn_Gly.mol2
Reading MOLECULE named GROup
> check LIG1
Checking 'LIG1'....

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Warning!
The unperturbed charge of the unit (-1.001000) is not zero.
Checking parameters for unit 'LIG1'.
Checking for bond parameters.
Checking for angle parameters.

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Error!
Could not find angle parameter: CA - CT - H1

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Error!
Could not find angle parameter: CA - CT - H1

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Error!
Could not find angle parameter: CA - CT - OS

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Error!
Could not find angle parameter: OS - C - N

/raid/data/pschrank/Anaconda3/anaconda/envs/gmxMMPBSA/bin/teLeap: Warning!
There are missing parameters.
check:  Warnings: 1
Unit is OK.

This ERRORS get only produced, when the topology for the mutant variant get generated but not for the wildtype (It seems it can't recognize Angle parameters for bonded atoms, that should correlate to the build in GLYCINE in Position 2?)

I don't if this helps you but I figured posting it doesn't make things worse :D

Best regards

Paul

[Paul Schrank]

Additionally I see in the leap input for the wildtype:

source leaprc.protein.ff14SB
source leaprc.gaff2
loadOff atomic_ions.lib
loadamberparams frcmod.ions1lm_126_tip3p
set default PBRadii mbondi2
REC1 = loadpdb _GMXMMPBSA_REC_F1.pdb
LIG1 = loadmol2 CBZ_Asn_Gly.mol2
loadamberparams _GMXMMPBSA_CBZ_Asn_Gly.frcmod
check LIG1
saveamberparm LIG1 LIG.prmtop _GMXMMPBSA_LIG.inpcrd
REC_OUT = combine { REC1 }
saveamberparm REC_OUT REC.prmtop _GMXMMPBSA_REC.inpcrd
COM_OUT = combine { REC1 LIG1 }
saveamberparm COM_OUT COM.prmtop _GMXMMPBSA_COM.inpcrd
quit

So the Amber parameters for the ligand get loaded in before the LIG check command. Whereas for the mutant:

source leaprc.protein.ff14SB
source leaprc.gaff2
loadOff atomic_ions.lib
loadamberparams frcmod.ions1lm_126_tip3p
set default PBRadii mbondi2
MREC0 = loadpdb _GMXMMPBSA_MUT_REC_F0.pdb
MREC_OUT = combine { MREC0 }
saveamberparm MREC_OUT MUT_REC.prmtop _GMXMMPBSA_MUT_REC.inpcrd
LIG1 = loadmol2 CBZ_Asn_Gly.mol2
check LIG1
loadamberparams _GMXMMPBSA_CBZ_Asn_Gly.frcmod
MCOM_OUT = combine { MREC0 LIG1 }
saveamberparm MCOM_OUT MUT_COM.prmtop _GMXMMPBSA_MUT_COM.inpcrd
quit

The Amber parameters for the ligand get loaded in after the check LIG1 command. Is this maybe related to the issue?

[Mario Sergio Valdes]

Yes, this is an inconsistency, but don't have any effect on the calculation because the fail is in the check command... Previously to saving the parametrized file, the frcmod is loaded, so everything works fine. In any case, I fixed and nothing changed. I trying with other systems to check if this happens. Do you have tested this system and mutation using the topology (-cp option) file instead of the structure?

[Paul Schrank]

Saldy I don't have... When I try running using the top file I always get the error: 

ParameterError: Not all improper torsion parameters found

I dont really know why this is happening, as I think I have all necessary Position restrains in my files (example for this system --> see in appendix

This is why I used the alternative method by providing the ligand parameters using the -lm option 

Do you know, what might be the issue?

[Paul Schrank]

Hello again! 

I managed to fix the error and the high delta H energies. After carefully analyzing my topology I found, that the error with the improper torsion parameters comes from the Histidines in my protein. Instead of HIP the double protonated Histidines were called HIS (which is no viable entry in the rtp of the amber14sb ff), so they were identified randomly as HID and HIE residues. So when trying to parse the topology, the improper torsions for these His residues couldn't be found as there is no entry to assign for protonated ND2 and NE1 using the HID or HIE entries. So if anyone else has issue related to this, please check the naming of your Histidines of your system and which naming describes which protonation state in your Forcefield! 

So before generating the topology I needed to rename the HIS to HIP - and voila the error vanished :)

Afterwards I run again the Glycine scan as kindly suggested by you using the cp flag and also voila the high delta H problem vanished! (you can check the FINAL_RESULT.dat if you want :D). So I guess the Ala scan might have some issues with calculating the internal bonding terms for the mutant receptor, when not using the top file (at least for my system). 

So thank you the issue got resolved!

Best regards 

Paul

[Mario Sergio Valdes]

This is good news.  I have been swamped and unable to review this case. Very sorry. I will try to solve it as soon as possible. The use of the structure to generate the Amber topology will be removed in future versions since it often poses this type of problem. The ideal would be to test this system with the topology. In any case, I'll try to cover all possible issues emulating this system.