result_scaled_water.txt (7.8 KB)
result_scaled_complex.txt (7.8 KB)
result_unscaled_water.txt (7.8 KB)
result_unscaled_complex.txt (7.7 KB)
I have attached the files.
I see that you used equilibrium FEP. Recently, there were some gromacs bugs fixed relating to the foreign lambda calculations, thus the results might be gromacs version dependent. I would recommend trying to run with the latest gmx version patch or try non-equilibrium switching.
I’m using Gromacs 2022.4
Okay, thanks.
And in Gromacs 2022.4, this has been addressed. https://manual.gromacs.org/2023-rc1/release-notes/2022/2022.4.html
Another thing is that this is a charge changing mutation which I am doing but I am keeping the charge of the system neutral by co-alchemical ion approach. Could this be be a problem?
Such a perturbation will be more difficult to converge and requires many lambda windows and longer sampling time. It seems that your simulations are very short. The standard deviation of the end windows is large. I would recommend to, firstly, calibrate your protocol against known protein-ligand benchmark set to ensure that you are able to reproduce those results.
1 Like
pmx ligandHybrid -i1 ./lig_1/lig1.pdb -i2 ./lig_2/lig2.pdb -itp1 ./lig_1/MOL.itp -itp2 ./lig_2/MOL.itp -pairs pairs2.dat
ligandHybridTop__log_> Reading ligand 1 from: "./lig_1/lig1.pdb"
ligandHybridTop__log_> Reading ligand 2 from: "./lig_2/lig2.pdb"
ligandHybridTop__log_> Reading topology 1 from: "./lig_1/MOL.itp"
ligandHybridTop__log_> Reading topology 2 from: "./lig_2/MOL.itp"
ligandHybridTop__log_> Reading file with atom pairs: "pairs2.dat"
ligandHybridTop__log_> Assigning forcefield parameters....
ligandHybridTop__log_> Superimposing mol2 on mol1
[13:07:42]
****
Range Error
atomId
Violation occurred on line 62 in file /home/conda/feedstock_root/build_artifacts/rdkit_1651075169244/work/Code/GraphMol/Conformer.cpp
Failed Expression: 30 < 30
****
ligandHybridTop__log_> ERROR: fitting mol2 on mol1 failed
ligandHybridTop__log_> Superimposing mol1 on mol2
[13:07:42]
****
Range Error
atomId
Violation occurred on line 62 in file /home/conda/feedstock_root/build_artifacts/rdkit_1651075169244/work/Code/GraphMol/Conformer.cpp
Failed Expression: 30 < 30
****
ligandHybridTop__log_> ERROR: fitting mol2 on mol1 failed
ligandHybridTop__log_> Initializing the main object for LigandHybridTopology
ligandHybridTop__log_> Starting hybrid topology generation
ligandHybridTop__log_> Making atom pairs.....
Traceback (most recent call last):
File "/home/ps004/anaconda3/envs/pmx/bin/pmx", line 8, in <module>
sys.exit(entry_point())
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/scripts/cli.py", line 97, in entry_point
PmxCli()
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/scripts/cli.py", line 44, in __init__
getattr(self, args.command)()
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/scripts/cli.py", line 64, in ligandHybrid
ligandHybrid.entry_point()
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/scripts/ligandHybrid.py", line 189, in entry_point
main(args)
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/scripts/ligandHybrid.py", line 308, in main
hybridTop._makeHybridTop( )
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/ligand_alchemy.py", line 2087, in _makeHybridTop
self._make_atom_pairs( )
File "/home/ps004/anaconda3/envs/pmx/lib/python3.10/site-packages/pmx/ligand_alchemy.py", line 2172, in _make_atom_pairs
a4 = self.m4.fetch_atoms(n2,how='byid')[0]
IndexError: list index out of range
And in another of the case where I have noticed ligA and ligB when mapped and hybrid is made it all works fine whereas if ligB and ligA are mapped and one tries to obtain the hybird there’s an error message with
ligandHybridTop__log_> Error: Something went wrong while assigning bonds! ligandHybridTop__log_> A-> Atom1: 5-C5' Atom2: 22-H02 ligandHybridTop__log_> B-> Atom1: 5-C5' Atom2: 21-H03
Regards,
Pallav
This means, you shouldn’t forget to provide ligandB as -i1 and ligandA as -i2 for ligandHybrid when using pairs2.dat
The error that you mention is, because of the same reason of using wrong inputs for pairs2.dat mapping
1 Like
Hi Vytas,
There’s something strange I have been noticing with regards to charge-changing perturbations.
I was intending to reproduce some data from this paper https://pubs.acs.org/doi/10.1021/acs.jctc.8b00825# and observed something strange with regards to atomMapping. As noted in this paper for charge changing perturbation from 39 to 34, I get almost similar ddG data as stated in the paper(about 0.4 Kcal/mol difference) but the problem is when I use atommapping and make hybrids from 34 to 39, I get somewhat different ddG(about 2Kcal/mol difference). Another thing I have observed with charge changing perturbation is that never dG’s of Ligand A to Ligand B and from Ligand B to Ligand A matches, however in case of charge conserving mutations, both these transformations are similar.
I would be very grateful if you could throw some light onto this.
Regards,
Pallav.
If so is the case, what is the difference between the pdb files mergedA and mergedB which I have noticed to give different structures?
As the “pmx -h” explains: mergedA is “hybrid structure based on the ligand 1”, merged B is “hybrid structure based on the ligand 2”, i.e. structure of ligand1 is used to build mergedA, while structure of ligand2 is used as a basis in mergedB.
Such a hysteresis in the results suggests a convergence issue. How large are uncertainties of the estimates? Are you using non-equilibrium free energy calculation method? If so, are the work distributions overlapping?
If you suspect a dependence on which ligand is used as a structural reference, maybe you are suffering from the lack of convergence in conformation conversion between ligand1 and ligand2?
result_39_34_water.txt (7.6 KB)
result_39_34_complex.txt (7.6 KB)
result_34_39_water.txt (7.6 KB)
result_34_39_complex.txt (7.6 KB)
I am using equilibrium FEP.
And that is particularly only when there’s charge changing perturbation. However for a charge-conserving perturbations no matter which way is the atom mapping done, the dG’s and ddG’s end up to be the same.
From the results it seems that 95% CI of the dG in complex covers range larger than 2 kcal/mol, so your results can only be interpreted up to this level of certainty. Charge changing modifications are more difficult to converge. I assume you are using a co-alchemical ion approach to keep the box neutral during the transition?
yes. co-alchemical ion approach.
Would the double box single system approach alleviate such problems?
Not likely, you are experiencing convergence issues. Sampling longer or using enhanced sampling could help. Alternatively, can instantiate end states with their respective conformers
1 Like
The paper quotes of about 10ns per lambda window. I find it strange that no one seeked the reverse transformations data.
They probably use HREX with REST to enhance sampling.
The problem is resolved.
For the AB hybrid, in the above mentioned case from 39 to 34, I was mutating a sodium ion to a dummy atom.
For the BA hybrid, all I had to do, was to grow a sodium ion from a dummy. I overlooked this aspect and was using a sodium and chloride ion respectively and turning the chloride ion to a dummy thereby the dG’s were not similar and strangely differing with large deviations.
Regards,
Pallav.