Cyclohexane parameterization

[ramess101]

@mrshirts @msoroush @jpotoff

As part of the JCED manuscript, we are demonstrating how MBAR can be used to parameterize cyclohexane (which currently does not have any Mie parameters). We performed simulations with the TraPPE force field and then predict VLE for a range of lambda, epsilon and sigma. With basis functions VLE properties can be predicted in less than a minute.

I decided to just perform a 2-D scan of epsilon-sigma for different values of lambda.

The plots below are not publication quality yet, e.g., the color scales only apply to the bottom right plot, the contours are not labeled, I only included a subset of the data to speed-up the analysis which caused the contours to be noisier than normal, and I need to scan higher value of epsilon for lam = 14-18. I will reanalyze higher epsilon today and work on cleaning up the plots, but I think they are good enough for discussion.

The top plot shows the scoring function (as suggested by Mick et al.) while the bottom plot is the average number of effective samples in the liquid phase.

The main takeaways are that we get pretty reasonable/smooth contours for lam = 12 even though TraPPE was our reference. We could try excluding liquid density and enthalpy of vaporization for state points where the liquid phase has fewer than 50-100 effective samples to remove any erroneous liquid phase MBAR estimates.

[ramess101]

@mrshirts @msoroush @jpotoff

We know that reweighting is problematic for predicting liquid densities when Neff is low (large changes in sigma and lambda). I tried removing all liquid properties with Neff < 50 from the scoring function, but this essentially changes your target properties, which leads to a very unsmooth surface. By contrast, we could ignore liquid properties even when Neff > 50 while we are exploring different lambdas and just compute the scoring function for vapor properties. In this case, we get a much smoother response surface, although the minimum is quite shallow.

[ramess101]

@mrshirts @msoroush @jpotoff

I repeated the scan over a higher epsilon region. The figures still need a lot of attention, but this is shows the region we will want to focus on:

[ramess101]

@mrshirts @msoroush @jpotoff

I have cleaned up this figure by processing all the snapshots. Also note that I decided to remove the stars that denoted the optimal parameter set, mainly because the uncertainties can impact what the true optimal is.

[msoroush]

@ramess101 it looks beautiful. Do you have phase diagram with optimum parameter for each lambda?

[ramess101]

@msoroush

I have not yet made such a plot, but I could easily do so.

Note that these results were obtained just by simulating the TraPPE force field. So we would want to resimulate the optimal eps/sig for each lambda to get a refined plot of the scoring function.

[jpotoff]

@ramess101 I'm getting jealous, Richard. You've taken our heat maps up a notch. These are really nice figures.

[ramess101]

Thanks @jpotoff You will notice that the version I sent you last night has a few outliers that I tracked down (my root solver still isn't perfect for equating the pressures sometimes). So the version here is a bit nicer.

[ramess101]

@msoroush

So plotting the phase diagram for the optimal eps/sig for each lambda does not look great. This is because MBAR does not predict the liquid densities accurately for lam > 12 (when the simulations are performed with TraPPE). The errors in rhol are essentially random, so all eps/sig are penalized equally and the optimization is basically just for the vapor phase.

Here are the percent deviations (top) and phase diagrams (bottom) for the optimal parameter sets and for TraPPE:

From this figure, the optimal eps/sig for lam =12 clearly cuts through the vapor phase data. Lam = 16 predicts rhov and Psat very accurately while the rhol and DeltaHv estimates are sporadic.

Again, we would want to simulate each of these optima to refine our search.

[jpotoff]

@ramess101 @msoroush Richard, do you want us to launch a grid of simulations around your optimal values so we can complete the parameter optimization?

[ramess101]

@jpotoff @msoroush

I don't think we necessarily need that for this manuscript (due in a week) since we aren't really reporting new MiPPE parameters. But we could certainly do that for a separate paper where we report optimal parameters for various cyclic compounds.

If we really want to complete the optimization, we could also just perform a single simulation for each lambda and then use MBAR to get the entire 2-d space like we already did for lambda = 12. But I don't think that we necessarily need to finish the optimization, since we already showed how good MBAR works for constant lambda and varying lambda.

However, it could be nice to compare the MBAR estimates with actually simulations for the parameter sets. How long would it take you to generate the 2-d grid if we were to do that?

[mrshirts]

If the manuscript is due in a week, then there really isn't that much opportunity to do a full parameterization; we can leave it as a demonstration of the approach.

@jpotoff, @msoroush, maybe now is the time to pause and say now that you've seen what can be done with this approach, what WOULD a full reparameterization paper look like? We don't need to overload everything in one paper.

[jpotoff]

@mrshirts @ramess101 Got it. Spooled up some simulations for the "optimum" lambda=16 for fun. Then found out Mohammad was already running a bunch of stuff...