GROMACS version: Any version
GROMACS modification: No
Dear everyone,
I am struggling to find some reference for equilibrium free energy simulations in case of annihilation (I mean: making vdWs and charges fade) of a, say, 50-atoms chromophore within a soluble protein cavity (i.e. in water), the protein being ~150 amino acids.
I understand that usually researchers are way more interested in finding hits or studying a lead phase for small molecules: 4_ns_production 10_ns_production
but this usually aims to qualitative trends and, moreover, considering small modifications, i.e. of part, of the original molecule.
Naturally, I wouldn’t dare to do a 4 ns (or the like) production per lambda step even if I were to adopt 20 lambda windows, in my case. Has anyone got any clue regarding at least simply the order of magnitude, like maybe 10-15 lambdas of 100 ns? Can’t really spot a reference… hope I have a “lead” or two :)
Thank you in advance!
P.S.: does nstdhdl = 10 , 100 or 200 make a difference? In this sense: say that you get the “true” free energy difference with 200, does this mean that lowering to 10 the keyword for printing values is lowering the accuracy of results (cause: autocorrelation) or, simply, it is useless and slowing down the simulation, but harmless?
In general, annihilating something as large as 50 atoms is challenging, especially if you want reliable results.
I won’t guarantee it will help, but I think you would have to aim at something in the range of 50 lambdas, 10 for Coulomb (might be possible to manage with less) and 40 for VdW (with tighter lambda distribution near the end points, especially the annihilated end point). You’ll definitely need soft-core interactions, but I would suggest using it only for VdW (you can explicitly set sc-alpha to 0 when you are only decoupling Coulomb). Regarding simulation time, I think you will have to experiment a little bit, but something in the range of 20-200ns per lambda point might be necessary.
If you are planning to analyse the results with BAR/MBAR check that you have good probability of sampling neighbouring lambda points. If there are problematic regions you can extend the simulation times of those.
Regarding nstdhdl, more data (samples) is good in general but, as you say, if the values are correlated it won’t help. I would say somewhere in the range of 20 to 100 is usually enough, and set nstcalcenergy to the same value.
You will not get lower accuracy of the results, but depending on how the uncertainties are calculated, the error might get underestimated if there is a strong autocorrelation.
thank you, once again. I think I am going to try even if that seems expensive, useless to say all your detailed information are precious for my current setting. Yes, I am using BAR, I rely on the printed relative entropies in these first trials to see at first glance if something I am doing is just way too coarse. Sincerely wish you a good continuation.
That sounds good. Just a final (for now) piece of advice: make sure that you are chaining your annihilation simulations, so that one starts from the end of the previous one. Otherwise, the equilibration time for the fully decoupled end point will be extremely long. You might already be doing this, of course.
Absolutely, you of course have many users to relate to but, hoping this also makes you even more passionate and satisfied about your didactic work in here, I remember our previous discussion and am sticking to the sequential runs accordingly: Solvation free energy calculation - #2 by MagnusL
Sounds good. I just realised that I, for some reason, read “50 atoms” as 50 heavy atoms. So my recommendation regarding number of lambda states (and probably simulation time per lambda state) was probably too high.