I’m trying to calculate the free energy of solvation of a particular molecule in water (amber03 and TIP3P parameters are used for solute and solvent, respectively). I’m ceasing electrostatic first and then van der Waals interactions in a protocol that includes: minimization, NVT (5ns), NpT (50ns) and Production (50ns). I have 31 lambdas and some calculations got a Fatal Error with this message:
“Fatal error:
There are 1 perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.2 nm, which is not supported. This can happen because the system is
unstable or because intra-molecular interactions at long distances are
excluded. If the latter is the case, you can try to increase nstlist or rlist
to avoid this.The error is likely triggered by the use of couple-intramol=no
and the maximal distance in the decoupled molecule exceeding rlist.”
I really do not understand this error. I checked some structures and did not see a problem that justifies, why some calculations finish without problems and others show this fatal error.
Did anyone face this problem? Should I extent my equilibration time before the production?
Thanks in advance,
Fred
This will fail as soon as some perturbed interactions exceeds this range. The rlist is normally automatically set by the Verlet-buffer-tolerance.
So as stated in the error, you can either extend your rlist or turn on the intra-molecular coupling.
I assume you don’t want to turn on any intra-mol coupling. So the solution would be to extend your rlist. Furthermore, to overwrite the auto rlist by Verlet-buffer-tolerance, you need to set Verlet-buffer-tolerance=-1. For me, the extension of my molecule reaches up to 1.6 nm, when it is fully extended. So I would recommend to check the size of your molecule.
Hi Tomo,
Thank you for your answer. It sounds like the best solution to overcome this problem. The average size of my molecules is smaller than yours.
Thank you again,
Fred
If we are doing drug like molecules, for example, ligpargen created an instance of ribociclib for me where the furthest atoms from each other are 1.754 nm apart. Should I be setting couple-intramol=yes or set rlist=rcoulomb=rvdw > furthest_points
As far I know, if you set your rcoulomb and rvdw (those are for the short-range electrostatic and short-range vdw) that large, your calculation will be dead slow (correct me if I am wrong Gromacs-experts). One has to be carefull with these settings.
I still see this with GROMACs version 2022.2 here https://ftp.gromacs.org/gromacs/gromacs-2022.2.tar.gz. The error returns exactly the furthest two points of the molecule you are trying to decouple. For instance if at step 100 you get this message:
Fatal error:
There are 1 perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.783 nm, which is not supported
or
Fatal error:
There are 1 perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.921 nm, which is not supported
This means, and you can confirm it with a simple distance algorithm using python itertools for all pairs of points, that the furthest two points of the molecule being decoupled would be 1.783 nm in the first case and 1.921 nm in the second case. Ironically nobody else tried to confirm this but me from what I have read.
I have tested for my system yesterday, and GROMACS 2022 still yielded to this error. So I believe that the assertion code is there and since your perturbed nb pair interaction exceeds your rlist, it threw this error.
I am having the same problem.
Even with the latest Gromacs, to calculate the free energy of a relatively large molecule (such as Fred’s molecule or bjwiley23’s Ribociclib) under couple-intramol=no conditions, I need to do rlist > furthest_points > rcoulomb=rvdw, which would require an enormous amount of time?
Does it mean that the calculation results are still not valid when I set rlist > furthest_points > rcoulomb=rvdw?
Also, if I set rlist > futhest_points, should the simultion_box_size > 2×rlist + futhest_points?
All requirements you list are really required. This is handled correctly in 2022.
But I would suggest to use couple-intramol=yes for large molecules, both to reduce the computational cost and for avoiding getting stuck in local minima when the molecule is decoupled.
Hi Hess,
Thank you for your answer. I am planning to perform TI on an elongated molecule called C12EO8 with a maximum molecular length of about 2.3 nm.
Is this molecule a “large molecule”?
I am also planning to calculate TI under the following conditions.
I understand that if I calculate with couple-intramol=yes, the Van der Waals and Coulomb interactions within the molecule will also be turned on and off, is this correct?
Thank you @hess for the response. Can you point to the relevant references to do the extra free-energy calculation, the implementation details and its importance?
Any serious text on binding free energy calculations should explain this. What you need to do is an identical free-energy coupling calculation but without solvent and subtract the resulting free-energy from the result from the calculation with solvent.
I encountered the same problem as @fredpontes20 with my ion channel-membrane system when running minimization step. Even after I set Verlet-buffer-tolerance = -1 and increased my rlist, I still get the same error despite how large rlist is. I also tried couple-intramol=yes and still get this error.
I’m not sure what’s happening here and what should I do to get this fix and get my minimization running.