ramess101 / MBAR_ITIC

Merging the MBAR and ITIC methodologies to optimize a TraPPE force switch model
0 stars 2 forks source link

Mie Optimization #4

Open ramess101 opened 7 years ago

ramess101 commented 7 years ago

Although this project was originally designed for the LJfs potential, I am also very interested in seeing how the Mie potential performs. Some disadvantages of using the Mie potential are:

  1. Gromacs requires use of tabulated potentials which take about 2-3 times longer than the native LJ
  2. Force switch is not allowed
  3. It is now a three dimensional parameterization with highly coupled parameters
  4. The lambda parameter should be an integer value to speed up computational costs (for Gromacs it is the same for integer and fractional lambda values because it is always tabulated, but for other software this might not be the case.)

The purpose of this issue is primarily to address disadvantages 3 and 4, since these directly impact the optimization approach. Specifically, this issue will demonstrate what we learn as we try to optimize a Mie potential for ethane using MBAR-ITIC.

ramess101 commented 7 years ago

One of our goals with MBAR (and PCF-PSO) is to rapidly convert the TraPPE (or other LJ based models) into a Mie based model. Here I show some of the difficulties in this endeavor.

First, the optimization itself is tricky. Leapfrog is supposed to scale well with number of dimensions but it also depends heavily on the bounds because if your initial players are spread out over a very large space it is difficult for them to converge. Furthermore, I believe leapfrog struggles with highly correlated parameters.

Second, MBAR is not great at predicting U and P when extrapolating from the TraPPE LJ model parameters to the Potoff Mie parameters.

In these figures we show the iterations used in the leapfrog algorithm in the 3-dimensional parameter space where light colored circles represent optimal points. The reference system is the TraPPE LJ model, which is 98 K, 0.375 nm lambda=12. The 'x' is difficult to see, not sure why it showed up as white.

image

The strange thing is that the optimization spent so much time in the upper left corner. This region does a decent job for rhoL but a very poor job for Psat. Although MBAR cannot find the exact optimum with a single reference system (i.e. the optimum region is quite different from the Potoff model), it does locate the general region where the optimum should be found. Meaning that it is likely that this will converge within a few more iterations.

It is worth noting that the SSE for rhoL of the optimum model is much worse than TraPPE and Potoff while the SSE for Psat is considerably better than TraPPE but slightly worse than Potoff. Again, we need to perform another iteration using this as our new reference system.

This figure shows what the SSE in U looks like (recall that the optimization currently only uses rhoL and Psat).

image

ramess101 commented 7 years ago

Currently I am repeating the optimization but using more informed bounds. Namely, epsilon can range from 98-130, lam from 14-18. Sigma will have the same bounds as before. Alternatively, I could let sigma be in a much smaller window and shifted towards higher sigmas. But I want to see if I can keep from including too much prior information.

ramess101 commented 7 years ago

Alternatively, we do not really need to start with such a poor guess. Instead, we could use the Potoff Mie potential and just improve it. Here I used the Potoff model as my reference and reoptimized the parameters to match rhoL and Psat, which were included in Potoff's optimization.

image

We do achieve a slight improvement but at the expense of using a lambda value very close to 17 (16.9). We could use an integer value of 17 to avoid fractional exponents, but still with combining rules a value of 17 could lead to fractional exponents. The point is, we can quickly optimize the Mie potential to improve an existing Mie model.

Here I show what U looks like although it was not included in the optimization. Clearly, the optimum for rhoL and Psat is slightly different from the optimum in U.

image

Again, we could improve the leapfrog approach by limiting the range of the parameters (especially lambda, since values less than 14 are known to be poor in this case).

ramess101 commented 7 years ago

In addition, we could keep lambda constant and just reoptimize epsilon and sigma. Here I have used Potoff's potential as the reference and reoptimized eps and sigma with lambda of 16. Notice that the optimum is very close to Potoff for rhoL and Psat. Again, U is not included in the optimization and demonstrates a slight disagreement with the rhoL and Psat optimum.

image

ramess101 commented 7 years ago

An advantage to using a fixed value of lambda and scanning the epsilon and sigma parameter space is that basis functions can be used to avoid needing to perform as many reruns.

ramess101 commented 7 years ago

Revisiting the issue of optimizing Mie parameters using TraPPE (UA-LJ) as the reference. Here we see an example where the leapfrog algorithm converged to a local minimum of sorts. It found a region with very good rhoL but very poor Psat.

image

image

ramess101 commented 7 years ago

I have repeated the conversion of the TraPPE-UA from a LJ potential to a Mie potential. I have increased the number of players in the leapfrog algorithm and performed five cycles where I talk the optimal epsilon, sigma, and lambda as the reference. Note that I forced lambda to be between 14-18 to try to expedite the optimization.

These leapfrog results use the TraPPE as the reference (Cycle 0):

image

image

These take the previous optimal as the reference (Cycle 1):

image

image

And so on...

Cycle 2:

image

image

Cycle 3:

image

image

Cycle 4:

image

image

One problem is that I believe I shrunk the boundaries too quickly. I thought this would speed up the optimization but I think I shrunk the bounds faster than the optimization was converging. For example, notice that cycle 2 converges in the corner.

That being said, I also think this is just a very difficult optimization to do in 3 dimensions with such correlated parameters. I think that it might be best to perform an individual optimization at each value of lambda instead. This would increase computational costs, but maybe we can use PCF-PSO to figure out which lambdas to try and to find the best reference parameters.

mrshirts commented 7 years ago

If there are correlated parameters, wouldn't it make sense to do some variant of gradient optimization? Was there an issue doing that?

On Thu, Aug 3, 2017 at 12:36 PM, Richard Messerly notifications@github.com wrote:

I have repeated the conversion of the TraPPE-UA from a LJ potential to a Mie potential. I have increased the number of players in the leapfrog algorithm and performed five cycles where I talk the optimal epsilon, sigma, and lambda as the reference. Note that I forced lambda to be between 14-18 to try to expedite the optimization.

These leapfrog results use the TraPPE as the reference (Cycle 0):

[image: image] https://user-images.githubusercontent.com/23408516/28936791-07b5b9be-7846-11e7-877a-bfb2cfbdae97.png

These take the previous optimal as the reference (Cycle 1):

[image: image] https://user-images.githubusercontent.com/23408516/28936804-117ac606-7846-11e7-89b3-93f0376a1336.png

And so on...

Cycle 2:

[image: image] https://user-images.githubusercontent.com/23408516/28936808-1720ada0-7846-11e7-8c7b-1c85e12616e4.png

Cycle 3:

[image: image] https://user-images.githubusercontent.com/23408516/28936828-266ed14c-7846-11e7-9ab0-0481d064eecb.png

Cycle 4:

[image: image] https://user-images.githubusercontent.com/23408516/28936832-2ad719e2-7846-11e7-9688-e5c39bcc8bbb.png

One problem is that I believe I shrunk the boundaries too quickly. I thought this would speed up the optimization but I think I shrunk the bounds faster than the optimization was converging. For example, notice that cycle 2 converges in the corner.

That being said, I also think this is just a very difficult optimization to do in 3 dimensions with such correlated parameters. I think that it might be best to perform an individual optimization at each value of lambda instead. This would increase computational costs, but maybe we can use PCF-PSO to figure out which lambdas to try and to find the best reference parameters.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ramess101/MBAR_ITIC/issues/4#issuecomment-320054374, or mute the thread https://github.com/notifications/unsubscribe-auth/AEE31AAgpurr8fhLr7lWPM7vWkM7Emg3ks5sUhNKgaJpZM4OkIlK .

ramess101 commented 7 years ago

If there are correlated parameters, wouldn't it make sense to do some variant of gradient optimization? Was there an issue doing that?

We have included the ability to use gradient based optimizers and they work very well if implemented properly. However, they are finicky to work with. The primary difficulty is knowing how much of change in epsilon, sigma, or lambda should be used to calculate the gradients. Since these gradients are obtained from simulation, the simulation noise may trump the variation caused by changing the parameters. In addition, most of the bounded scipy optimizers require evaluations of the objective function at the boundaries, where MBAR will be less accurate (assuming the reference is in the center of the bounded parameter space). So, in brief, we could use the gradient based approaches but they require a lot of hand-holding. Also, gradient approaches actually converge to the same result as the leapfrog approach. The reason I like leapfrogging is that it still samples from a wide parameter space to truly demonstrate that we have found the global minimum.

So I think it is less of an optimization issue as it is an issue of how quickly MBAR will transition from lambda = 12 to lambda = 16. It just seems to not want to make a big jump. For example, starting with LJ the new optimum (regardless of how it is obtained) will have a lambda around 12.1, which shows how long it will take for lambda to get to 16 (if it will ever get there). For this reason, it might be ideal to just have a few independent optimizations with different values of lambda where MBAR just optimizes epsilon and sigma for a constant lambda value. In other words, find the optimum for lambda of 12, 14, 16, 18 and then compare the objectives. This would also facilitate using basis functions as lambda would be held constant from one reference simulation to the MBAR reweights.

It is worth pointing out that we are using rhoL and Psat as our targets from the get go. The initial reference (TraPPE) gives a very good rhoL with a poor Psat. rhoL is related to Z vs 1/T (on the isochore) which means that it all comes down to P from MBAR. P is not predicted as accurately as U for different parameter sets from the reference. So, perhaps for early cycles we could use U to locate the optimal region and then transition to rhoL and Psat. But if we are going to do this we might as well just use PCF-PSO since we know that works very quickly with U.

ramess101 commented 7 years ago

One thing that Richard and I noticed is that MBAR appears to be more accurate at predicting U than P for large differences in the parameters, specifically sigma and lambda. So I decided to optimize just for U first when making the large jump from LJ to Mie and it worked a lot better. Here is the first rerun cycle:

image

Notice that epsilon is about 121.7 (where Potoff has 121.25) and lambda is around 15 while sigma stays fairly constant. Here is the second rerun cycle:

image

Notice that the epsilon is still in the ballpark (like 119) while lambda continues to increase (it is 15.5) and sigma has started to take a hit (around 0.381) but it does need to increase to around 0.3783 nm. The sigma of 0.381 nm will lead to poor P on the isochores (which will cause poor rhoL, Psat) but currently we just want to get U good before we try to do both U and P.

Here is the third rerun cycle:

image

Notice that epsilon is still around 120, sigma is up to 0.384 and lambda to 17.3. This demonstrates that it might be useful to include a constraint on sigma and rmin.

Here is the fourth rerun cycle:

image

Here is the fifth rerun cycle:

image

Notice that there is a wide range of sigma values that match U fairly accurately, (but some of these will likely have very poor pressures). Because of the low precision in sigma with U it is difficult to use a convergence criterion that is based on the values of the parameters from one rerun cycle to the next. Instead, it may be better to use a convergence criterion based on the objective function. That is, once the objective function is not improving to within a given tolerance the system is said to have converged. Then, we can do U and P simultaneously (since MBAR's prediction of P should be accurate in that region). Then, we could transition to the more traditional rhol and Pv optimization.

Again, a simpler alternative is to just optimize eps and sigma for a given lambda. Originally I thought this would require more simulations, like a reference simulation for each lambda. But if we start with U we could just have lambda be an integer value in the leapfrog algorithm. By using U we should still get to the same optimal lambda and then we do not need to worry about having a different reference at each lambda.

ramess101 commented 7 years ago

If there are correlated parameters, wouldn't it make sense to do some variant of gradient optimization? Was there an issue doing that?

One more point regarding the use of leapfrog algorithm. Leapfrog is very nice for using constraints that a given parameter is integer. I was unsuccessful at implementing a constraint for integer parameters using any of the scipy optimizers. Perhaps Michael knows how to do this, but when I searched online I found a few discussions stating that there are no good scipy optimizers to handle this. By contrast, with leapfrog you simply require that every leap of lambda remain integer.

But this begs the question, for Mie potentials, how important is it to use integer exponents? The main benefit is computational costs. Since Gromacs does not have Mie potentials I am using tabulated potentials already. So the computational cost is the same on my end regardless of if it is an integer exponent. But for users of the force field I have noticed a considerable decrease in speed when using floating point exponents. That being said, although the Gromacs reruns will cost the same, if we use a constant (integer) value for lambda we should be able to use basis functions to speed up MBAR, right Michael? Which might be a really big benefit. Do you (Richard and Michael) have an opinion regarding whether we want to constrain the optimization to use integer exponents? I am leaning towards yes, I just want to see if we are in agreement.

Also, in a related topic, we might find that the Exp-6 model is better than the Mie potential. I am anticipating this to be true at higher pressures. In that case, the Exp-6 model uses a parameter (typically referred to as alpha) that is analogous to lambda in the Mie potential, i.e. it controls the repulsion. The Exp-6 is already much more computationally demanding. Do either of you know if it matters if alpha is an integer for the Exp-6? I believe it does not matter since alpha is multiplied by the distance which is inside the EXP. So I think the EXP(alpha*(1-r/rmin)) will have the same cost regardless. But I believe the Exp-6 model cannot be used with the basis function approach (to speed-up MBAR), right Michael?

mrshirts commented 7 years ago

The Exp-6 is already much more computationally demanding Though it depends on the context; if tabulated, not very much. And I thought that exponential functions are relatively well implemented in most languages.

Do either of you know if it matters if alpha is an integer for the Exp-6 Alpha is generally a free parameter.

My instinct for Mie is to treat lambda as a free parameter for optimization, and then one can backtrack afterwards and use whatever integer value is closest. But it really depends on what the scientific goals are.

In terms of basis function for reweighting; with both Buckingham and Mie, you can handle two of the parameters with basis functions, but not the third. For Mie, you have essentially two basis function associated with 1/r^n and 1/r^m, for Buckingham one with exp(-Br) and one with 1/r^6. So you have to reevaluate direct energies with each choice of lambda (n,m) or B, but then the other 2D can be filled in with basis functions.

ramess101 commented 7 years ago

Though it depends on the context; if tabulated, not very much. And I thought that exponential functions are relatively well implemented in most languages.

True, most languages have the exponential function well implemented, but it is still more expensive than the LJ or Mie if they are native to the simulator. But like you said, if the Mie is tabulated (like in Gromacs) the Exp-6 might actually be about the same.

ramess101 commented 7 years ago

In terms of basis function for reweighting; with both Buckingham and Mie, you can handle two of the parameters with basis functions, but not the third. For Mie, you have essentially two basis function associated with 1/r^n and 1/r^m, for Buckingham one with exp(-Br) and one with 1/r^6. So you have to reevaluate direct energies with each choice of lambda (n,m) or B, but then the other 2D can be filled in with basis functions.

OK, let me clarify my question regarding the Mie potential. First, we are going to assume m=6 (this is pretty common for optimizing Mie potentials, although we could look at relaxing this if needed). So the 1/r^6 term can work with basis functions. I was thinking that if we limit lambda (or "n") to be integer values then we could use basis functions for a given lambda. For example, we use the basis function approach to obtain the entire epsilon and sigma parameter space when lambda = 12, then we have a separate basis functions when lambda = 14, lambda = 16, etc. This would only be possible using integer values of lambda. Does that make sense? Is that a feasible approach?

Let me make sure that I understand the basis function approach. My understanding is that you just need to perform a few reruns (with no repulsion/attraction or something like that to get the array of 1/r^n and 1/r^m) and then you can calculate the energies for any epsilon and sigma (assuming lambda, or n and m, are constant). Is that correct? By not needing to perform a rerun at every value of lambda, which is necessary when lambda is not integer, we could speed-up the MBAR process considerably. Currently it takes about 30 seconds to do a rerun at the 19 ITIC state points (for a new epsilon, sigma, lambda). Basis functions would be more like microseconds, right? So if we developed basis functions for each value of lambda we could probably perform more rigorous optimizations or Bayesian methods.

Also, this would allow us to investigate anisotropic (AUA) models. In other words, we can repeat this process for different bond lengths. This may be useful for ethane and other diatomics. But for larger compounds we are unlikely to use AUA.

ramess101 commented 7 years ago

I was also thinking that it could be beneficial to constrain the optimization when only using U as the target property. Since sigma (and rmin) do not affect U as much as P, we don't want sigma to become non-physical by only optimizing to match U. This was a common problem for PCF-PSO when only using U. The solution I found is to use a constrained optimization for the early runs (i.e. when epsilon and lambda are changing considerably). The constraint is just that if lambda > 12 then sigma > sigma_TraPPE and rmin < rmin_TraPPE and the inverse if lambda < 12. This plot shows what this constraint looks like for lambda >= 12:

image

Each line (or point in the case of LJ lambda = 12) corresponds to the sigma values that are acceptable for a given lambda value. Notice that they are within the constraint that sigma > 0.375 and rmin < 0.4209.

With this constraint, higher values of lambda have a larger range of possible values of sigma. This can lead to a non-physical bias that favors higher lambdas and can prohibit the optimizer from finding the true global minimum. However, the leapfrog algorithm can apply this constraint such that the probability of obtaining a given lambda is constant. This is done by only changing sigma if an attempted move violates the constraint while leaving lambda as the proposed value. Again, this is just to help the optimizer find the true global minimum.

mrshirts commented 7 years ago

Let me make sure that I understand the basis function approach. My understanding is that you just need to perform a few reruns (with no repulsion/attraction or something like that to get the array of 1/r^n and 1/r^m) and then you can calculate the energies for any epsilon and sigma (assuming lambda, or n and m, are constant). Is that correct? By not needing to perform a rerun at every value of lambda, which is necessary when lambda is not integer, we could speed-up the MBAR process considerably. Currently it takes about 30 seconds to do a rerun at the 19 ITIC state points (for a new epsilon, sigma, lambda). Basis functions would be more like microseconds, right? So if we developed basis functions for each value of lambda we could probably perform more rigorous optimizations or Bayesian methods.

Correct. My point is that even if you have to reevaluate lambda, then you can use basis functions for the other dimensions, which should still save time.

In other words, we can repeat this process for different bond lengths.

And I think mapping could be useful here.

ramess101 commented 7 years ago

Correct. My point is that even if you have to reevaluate lambda, then you can use basis functions for the other dimensions, which should still save time.

Yes, I think this would be very useful. But recall that this requires using integer (or at least specified) lambda values.

And I think mapping could be useful here.

Yes, we have yet to explore how well we can do this.

ramess101 commented 7 years ago

Here I have limited lambda to integer values. Admittedly this plot is not the most useful. It would be better to plot epsilon vs sigma for each value of lambda. The optimal (the light colored region) is in the lambda = 14 plane. This is interesting since the optimal with non-integer lambdas was about 14.9. I am a little concerned about this. Perhaps the optimization did not work properly.

image

ramess101 commented 7 years ago

Here I am using integer values of lambda with a constraint on sigma that requires sigma > sigma_TraPPE and rmin < rmin_TraPPE when lambda > 12 (which is always the case here). Here is how the leapfrog algorithm results look:

image

Again the optimum has lambda of 14. Epsilon is pretty good (around 116). And since sigma is constrained it stays in a good region as well. Though, actually, if you notice the axis range for sigma in this plot is shifted compared to the plot in the previous comment. That plot used unconstrained optimization. I guess I should have just checked to see if the sig_opt was in the constrained region. So that means that for the first rerun cycle the constraint on sigma doesn't matter much. But the higher rerun cycles were when we saw sigma go crazy.

To verify that this is actually finding the global minimum (when lambda is an integer) I am going to optimize lambda for lambda = 13, 14, 15, 16, 17, 18, etc. and see if lambda = 14 is truly the optimal.

mrshirts commented 7 years ago

Yes, it does require specified values. No need to be integer.

On Tue, Aug 8, 2017 at 8:35 AM, Richard Messerly notifications@github.com wrote:

Correct. My point is that even if you have to reevaluate lambda, then you can use basis functions for the other dimensions, which should still save time.

Yes, I think this would be very useful. But recall that this requires using integer (or at least specified) lambda values.

And I think mapping could be useful here.

Yes, we have yet to explore how well we can do this.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/ramess101/MBAR_ITIC/issues/4#issuecomment-320975352, or mute the thread https://github.com/notifications/unsubscribe-auth/AEE31LGkahaBd147wXnX1ncCyz-MpdTCks5sWHI4gaJpZM4OkIlK .

mrshirts commented 7 years ago

So, I think generally, if it's matter of answering the problem for 2-3 systems (or even 5-10), then I don't think we need to spend too much time trying to find the best approach -- iterative optimization seems to be working fine to narrow things down. If you want to use this over lots of different systems, then we should back up and strategize optimal approaches from the beginning (back up, enumerate all the things that have work and haven't, and stitch back together a full plan).

On Tue, Aug 8, 2017 at 2:04 PM, Michael Shirts mrshirts@gmail.com wrote:

Yes, it does require specified values. No need to be integer.

On Tue, Aug 8, 2017 at 8:35 AM, Richard Messerly <notifications@github.com

wrote:

Correct. My point is that even if you have to reevaluate lambda, then you can use basis functions for the other dimensions, which should still save time.

Yes, I think this would be very useful. But recall that this requires using integer (or at least specified) lambda values.

And I think mapping could be useful here.

Yes, we have yet to explore how well we can do this.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/ramess101/MBAR_ITIC/issues/4#issuecomment-320975352, or mute the thread https://github.com/notifications/unsubscribe-auth/AEE31LGkahaBd147wXnX1ncCyz-MpdTCks5sWHI4gaJpZM4OkIlK .

ramess101 commented 7 years ago

In the near future we are mainly concerned with n-alkanes. This represents about 12 systems since we are only considering compounds that are found in REFPROP. That being said, I would say that we do want to have this work for a lot of different systems. Originally our focus is non-polar compounds, so most likely the lessons we learn for n-alkanes will transfer over.

I agree that we should back-up and strategize some more before pressing forward. Would you want to meet sometime to discuss this?

ramess101 commented 7 years ago

Here I verify that the true optimal (when reference is TraPPE and integer values of lambda) is at 14, although 15 has a very similar optimum. I admit that the following plots are not as helpful as they could be. I didn't think it merited my time to combine the SSE plots all onto one figure and have clear labels, etc. The main point is just that I verified the leapfrog with integer lambda and constraints converged to the right result.

Here I show the SSE as a function of leapfrog iterations when lambda can be any integer (it converges to 14)

image

Here I show the same thing but where lambda is forced to be 13:

image

Notice that the optimal is much higher than for lambda of 14.

Here I show the optima from 15-18 (i.e. the first ~75 iterations are 15 the next set are 16, etc.):

image

Again, notice that the optima have a higher objective function. Lambda of 15 has a very similar optima to lambda of 14 though.

Here I show the parameter space leapfrog results. This plot is for lambda 15-18:

image

This is for 13:

image

ramess101 commented 7 years ago

I decided to see what happens if I start with the Potoff model as the reference system. Interestingly, if I just use internal energy as the target it converges within 2 iterations. The optimal lambda is 17 although a lambda of 16 is very similar. It goes from 16 to 18 to 17 and stays at 17.

These three plots represent the rerun cycles. So the first reference system is Potoff and we scan lambda of 13-18 with constraints on sigma. Then we take the optimum as the new reference, etc. The lambda values increase from left to right so whenever SSE goes back up it is a new lambda value starting the leapfrog optimization.

image

Notice in this first rerun that the middle group (for lambda of 16) is slightly higher. I believe this is because the reference system has lambda of 16 so these MBAR predictions are actually more accurate. Whereas the other lambda values can find an eps/sig that has a lower U because they essentially just find a set with poor overlap and that allows for a single snapshot to dominate and give a better U.

image

image

Notice again the lambda of 17 is the final value but it is very similar to lambda of 16.

This is what the objective function looks like in 3-dimensions for the 3 rerun cycles. The strong frontier with respect to sigma and lambda is due to the constraint.

image

image

image

Here I have taken cross-sections for each value of lambda (sorry that the sigma axis is messed up, I need to fix that):

image

image

image

Recall that we only use U for the objective because P has more bias when predicted with MBAR. This is what the SSE P looks like with the sampled states:

image

image

image

Notice that in all three cases the optimum set for P is at lower lambda (at least when we consider the regions that were sampled for optimizing U). It is likely there is some eps and sig that would do better for P with each lambda.

P is directly related to rhol. Psat is a combination of the knowledge from U and P. So looking at Psat we can see how well MBAR would work with a single reference:

image

image

image

None of these sets seem to predict Psat very well. So we would clearly want to include Psat in our objective eventually.

ramess101 commented 7 years ago

I just realized that in all of these plots the epsilon axis is wrong. It should say epsilon (K) not (kJ/mol). Hopefully that was already obvious. It is not practical to go back and fix this now.

ramess101 commented 7 years ago

I have been testing the PCFR-ITIC approach (PCFR is pair correlation function rescaling). Currently, the best version of PCFR is PCFR-PMF where PMF denotes that we are rescaling the PCF based on the PMF. Namely, the PMF of the reference system is shifted by delta_u (the difference in pair-wise energies between reference and target systems). Although I have observed some deficiencies in this approach, specifically with large changes in sigma, it actually does a pretty good job at extrapolating to different lambdas (which might be most important). By this I mean that it provides reasonable estimates for rhol and Psat at lambda of 16 when the reference is the TraPPE LJ model. Here are the PCFR-ITIC contours for rhol.

First, for LJ (lambda = 12):

image

Now for lambda = 16:

image

Also, here is Psat for lambda = 16:

image

Although the trends are not perfect, notice that an optimization of rhol and Psat would likely result in parameters very close to Potoff's (121.25 and 0.3783). The sigma value is likely to be too low, but considering how you could get a decent estimate in one iteration is very encouraging.

By comparison, MBAR-ITIC struggles at predicting rhol and Psat when going from the TraPPE model to a lambda =16 model. This is primarily due to the low number of effective samples (poor overlap). Here is the average number of effective samples for the LJ model where the reference is the middle of the plot:

image

To clarify, yellow denotes 1000 effective samples whereas purple denotes ~1 effective samples.

Here is for lambda 16:

image

Notice that the large overlap is at lower sigma values. This is undesirable because we know that the optimal sigma should be at slightly higher sigmas for larger lambdas. Namely, for Potoff sigma is 0.3783 at the lambda = 16. This plot shows that Potoff (121.25, 0.3783) has ~1 effective samples:

image