Sampling issue for my system

[Jay-0520]

Hi there, I know sampling is not easy for FEP, and it seems I got a sampling issue for my systems. I'm trying to run a3fe for 100 compounds against 1 target protein. Can I ask two questions?

1. for one of my lambdas, I got the following stats

    ================================================================================
    FINAL SUMMARY: Predicted Maximum Efficiency Runtime vs Time
    ================================================================================
    Time/Rep(ns) Total(ns)  Actual(ns)  Predicted(ns)  Inter-run_SEM  Norm SEM      
    --------------------------------------------------------------------------------
      0.2          1.0        1.000             3.692          1.889          0.152046    
      0.4          2.0        2.000             5.721          2.069          0.235568    
      0.8          4.0        4.000             8.187          2.094          0.337114    
      1.6          8.0        8.000             11.407         2.063          0.469720 

when I extend the simulation length, the normalized SEM and predicted runtime also increase. As the result, it takes a long time to complete simulations for this lambda (has to be run for the maximum time, 30ns per simulation). Can I ask what is the recommended way to handle this case? would it be helpful if I create more replicas per lambda? I noticed that for the 5 repeats of this lambda, the corresponding statistical inefficiency values are 1.51, 3.27, 51.57, 3.97, and 6.54 at 1.6 ns.

2. I noticed that statistical inefficiency is not used to correct sems_inter_delta_g in process_grads.py. And by default, normalized SEM is computed by calling gradient_data.get_time_normalised_sems(origin="inter_delta_g"). Can I double check if this is designed to be this case? Besides, can I use other SEM metrics to handle the issue in question 1?

Thank you very much.

Best,
Jay Huang

[Roy-Haolin-Du]

Hi Jay, thanks for your question, it’s a really good and we can explore and learn more about it together!

1. For the question about normalized SEM and predicted runtime grow when you extend runs, based on my previous experience, this behavior is expected: it tends to increase initially, but later it stabilizes and may even decrease slightly.
   SEM × sqrt(T × R) makes the normalized SEM grow ~sqrt T , so the predicted optimal runtime increases firstly. But, later once the replicas start sampling equilibrium distribution Var_r(g̅r) ~ 1/(RT)，hence SEM = sqrt(C/(RT×R)) × sqrt(TR) = sqrt(C/R),is a constant with the central limit theorem satisfied.

   a3fe/a3fe/analyse/process_grads.py

   Lines 117 to 125 in 4e25894

   	squared_sem_inter = _np.var(means_intra) / len(run_nos)
   	mean_overall = _np.mean(means_intra)
   	# Store the final results, converting to arrays for consistency.
   	tot_sem = _np.sqrt(
   	squared_sem_inter + squared_sem_intra
   	) # This isn't really a meaningful quantity
   	sem_intra = _np.sqrt(squared_sem_intra)
   	sem_inter = _np.sqrt(squared_sem_inter)

   
   Shown one example I ran adaptive, the inter-run fluctuations are larger at the beginning and the differences increase, but at later stages they tend to stabilize and may even converge. Usually, it won’t go up to 30 ns per simulation.

2. The reason we use inter_delta_g is for the final ΔG directly. More care about differences between replicas impact on the final ΔG and concerned about the correlation between repeats. Therefore, allocating more time to λ’s with large inter-replica disagreement. Compared to the autocorrelation analysis of a single trajectory, there are some discussions on subsampling that you can take a look at. Should Subsampling be Recommended? choderalab/pymbar#545.

My recommendation: if you have successfully run other compounds, then continue running and see whether the SEM remains unchanged.
Please let me know if you have any thoughts. And Hi Finlay@fjclark, if you have time, I’d be very interested to hear your expert perspective as well. Many thanks.

[Jay-0520]

Thank you very much for the information and the detailed explanation @Roy-Haolin-Du.

I'm extending the simulations for now and I will come back with some analysis and stats (hope they can be useful). I also have some other questions related to sampling, e.g., the initial ligand is very mobile. But I will dig a bit more before bugging you here.

Really appreciate your help.

Best,
Jay Huang

[fjclark]

Hi both, thanks for the clear question and thanks for addressing this @Roy-Haolin-Du. What kind of overall errors are you getting @Jay-0520? To add a couple of comments:

when I extend the simulation length, the normalized SEM and predicted runtime also increase.

The discussion around table 2 here might be helpful. As Roy says, if the central limit theorem were satisfied, SEM would decrease with 1/sqrt(t) and hence the normalised SEM (SEM * sqrt(t)) would be a constant. This usually seems to be the case for the free leg, but is almost never the case for the bound leg due to bigger sampling issues (again see table 2). This is unfortunate, because the sampling time allocation algorithm is derived assuming that the central limit theorem holds. As you say, if each run is stuck in a local minimum, the inter-run SEM will stay the same, and the normalised SEM and predicted optimal runtime will both increase. There are three regimes:

a) Timescale of sampling issue << less than typical runtime: Algorithm works as expected, but results would still be good without it, so not really needed
b) Timescale of sampling issue ~ typical runtime: Sweet zone where the algorithm can help substantially compared to non-adaptive runs
c) Timescale of sampling issue >> typical runtime: Annoyingly common (pretty much all bound legs) and problematic zone which you're in, where the algorithm doesn't really help because it can't allocate enough sampling time to help. The only solutions are to run orders of magnitude more sampling, or to use better sampling (e.g. REST2, GCMC for waters, etc.).

Because we often encounter c), the adaptive algorithm only occasionally helps, as discussed in the paper. Also as discussed with table 2, in case c) the allocated run time will depend on lots of things which it should theoretically be independent of, such as the number of replica runs.

Can I ask what is the recommended way to handle this case?

I'm afraid I don't have any great solutions, but I would suggest

a) Understand your sampling issue. Isolate the lambda windows which are contributing most to the error over the whole calculation (e.g. look at the PMFs and see where they diverge, look at the Gelman-Rubin diagnostic for each window, look at the summaries of SEM as you're already doing) and do some analysis of these. For example, watching the trajectories is very often informative, and you can also look at the RMSDs and water in the binding site e.g. with calc.analyse(plot_rmsds=True) and this function might be handy. Also note that if you're running charged ligands you'll likely encounter much larger errors.
b) See if you can figure out how to address the sampling issue. This is tricky, because it's often ligand-specific and a3fe doesn't currently support replica exchange/ GCMC/ etc. Sometimes if you have slow changes between binding poses you can run a separate simulation for each pose (where the restraints isolate each pose) and combine the results. http://pubs.aip.org/aip/jcp/article-abstract/125/8/084901/916851/Parallelized-over-parts-computation-of-absolute?redirectedFrom=fulltext . You can also try more replicates -- the Coveney group recommends that for a given simulation time, you use the maximum amount of replicas where the simulation time remains long enough for proper equilibration. How useful this recommendation is in practice I don't know, because it's impossible to know the equilibration behaviour ahead of time. However, if you seem to have equilibrated (no overall slope to your convergence plots) with your current simulation time, you could decide on a reasonable fixed simulation time and rerun this with more repeats. This also follows the idea of https://onlinelibrary.wiley.com/doi/abs/10.1002/pro.5560070314

normalized SEM is computed by calling gradient_data.get_time_normalised_sems(origin="inter_delta_g"). Can I double check if this is designed to be this case?

This is intended. You can calculate the SEM both by looking at the deviations between repeats (effectively assuming each run is uncorrelated, but not worrying if the samples with each run are correlated), or by estimating the number of uncorrelated samples within each run (e.g. estimating the statistical inefficiency) and using this to calculate the SEM from the within-run data. These estimates should both agree, but where there are sampling issues, inter-run SEM > intra-run SEM. This is the basis of the Gelman-Rubin diagnostic in MCMC. We use the inter-run SEM because it tends to be less of an under-estimate, but note that both are usually underestimates and the choice of starting point e.g. input structure (even with ensemble equilibration) means that the repeat run results will still be correlated to some extent.

Besides, can I use other SEM metrics to handle the issue in question 1?

Unfortunately, I don't think using other metrics will help -- the inter-run SEM is the most realistic estimate of the SEM (even though still an underestimate). Using the intra-run SEM will give you less allocated sampling time due to being more of an underestimate, but this is problematic given more sampling is obviously required. There's not really an easy way round using a) better sampling (which is unfortunately not available in a3fe, although it's hard to know what exactly would improve your results wihtout understanding the sampling issue) b) orders of magnitude more sampling, which is likely unaffordable. Using more replicates might be a reasonable approach though.

Hope that makes sense/ is useful. Good luck with your system -- ABFE is hard...

Best,
Finlay

[Jay-0520]

Thank you @fjclark and @Roy-Haolin-Du. Sorry for my late response. I'm still trying to gather results. I have two questions about 1. resuming the calculation from where it stopped, and 2. running SOMD simulations concurrently in one GPU. Can I ask these questions here or I should open another ticket? Thanks.

[Roy-Haolin-Du]

Hi @Jay-0520,

Thank you for continuously exploring the applications of A3FE. We’re really glad to work together on solving these problems.
Please feel free to post your question here firstly. So we can have a quick discussion. If it turns out to be a more general issue, we can then escalate and broaden the discussion individually.

Thanks again.

[Jay-0520]

Thanks @Roy-Haolin-Du. The broader context is that I'm trying to run a3fe in one node (1 GPU + few GPUs) by using this script. When I looked at the codebase, it occurs to me that by default a3fe relies on slurm jobs for SOMD runs. What is your recommendation on running SOMD concurrently on one GPU and running MBAR in parallel on multiple CPUs? My script is about ~2000 lines (it works now), although most of them are about patching functions in a3fe, and I was wondering if there are better ways.

I found it less straightforward to resume the calculation from where it stopped. When my calculation get killed, it often results in inconsistent runtime for different runs (e.g., run_nos, 1, 2, 3, 4, 5,). It causes an issue for efficiency and equilibration checks because they need runtime to be consistent. So I need to fix runtime inconsistency firstly prior to resuming. What is your recommended way on fixing this? My approach is simply deleting the runs of a lambda that have different runtimes and rerun them from beginning.

Thank you very much.

[Roy-Haolin-Du]

I found it less straightforward to resume the calculation from where it stopped. When my calculation get killed, it often results in inconsistent runtime for different runs (e.g., run_nos, 1, 2, 3, 4, 5,). It causes an issue for efficiency and equilibration checks because they need runtime to be consistent. So I need to fix runtime inconsistency firstly prior to resuming. What is your recommended way on fixing this? My approach is simply deleting the runs of a lambda that have different runtimes and rerun them from beginning.

Hi, for resuming calculations, you are right, the key thing is that A3FE expects all repeats of all simulations to have the same total simulation time. To identify any inconsistent trajectories. You can use a quick check like:

for leg in calc.legs:
    for stage in leg.stages:
        for lam_win in stage.lam_windows:
            times = [sim.tot_simtime for sim in lam_win.sims]
            if not all(t == times[0] for t in times):
                print(lam_win.base_dir)
                print(times)

If you’re running adaptive, make sure all simulations to have the same total simulation time before continuing by using:

for leg in calc.legs:  
    for stage in leg.stages:  
        for lam_win in stage.lam_windows:  
            sim_times = [sim.tot_sim_time for sim in lam_win.sims]  
            max_sim_time = max(sim_times)  
            for sim in lam_win.sims:  
                if sim.tot_sim_time < max_sim_time:  
                    sim.run(simtime=max_sim_time - sim.tot_sim_time)

After that, you should be able to continue with:
calc.run(adaptive=True, runtime_constant=0.0005)

But, if your trajectory times look inconsistent (e.g. 200, 400, 600, 1200 … missing 800, 1000), or if a job was killed mid-step leading to “weird” run times, the safest fix is to wipe the job and restart cleanly:
sim.clean() then run a desired time.

Hope it helps!

[Jay-0520]

Yes, agree. I think in some cases, I have to delete the s3 file before resuming because resumed run will always start from this checkpoint, and resulting in a "gap" in simfile.dat. I agree the safest way is to re-start the job from the scratch.

[Jay-0520]

Okay, I will continue to run my simulations (my script seems to work now) and bring back some analysis later for this issue. Thanks!