To demonstrate why the Mie 16-6 model over-estimates viscosity at high densities at 293 K but not at the triple point, we are investigating the RDFs of propane at these two state points:

Clearly, the high pressure simulation has much closer CH3 and CH2 interactions. We are working on making these look a little prettier. @Kaiveria can you make the CH3 x CH2 RDF as well?

@mostafa-razavi

I think it would be valuable to test how sensitive our viscosities are to the bond constraints. We can use the harmonic bond potential recommend by de Pablo and compare with what we have for LINCS order = 8.

@jpotoff @msoroush @mostafa-razavi @jrelliottoh

In a previous study we developed joint distributions for the eps_CH3, sig_CH3 and eps_CH2, sig_CH2 parameter sets. These various combinations are denoted as MCMC parameter sets because we used Markov Chain Monte Carlo to sample from the posterior distribution using Bayesian inference. The results from this analysis are shown below (note that I have cut out parts of these figures that would only distract from our discussion):

The question remains as to how we want to quantify the uncertainty in CH and C. Although I could perform this Bayesian analysis, this would take some considerable effort and might not be the best approach anyways (different objective function, data sets, compounds, etc.) Rather, I propose that we just use the uncertainties as they have been provided in Mick et al. Here are the heat maps for these site types:

From these heat maps I can construct an approximate multivariate normal distribution. All I need to do is assign the variances in eps_CH, sig_CH, eps_C, and sig_C, along with their corresponding covariances. This is actually quite easy, although it requires some decision making on our part. This type of uncertainty is somewhat of a Type A and Type B because we are using statistical measures (the scoring function) but we are not using statistical methods to determine which scoring function values we would consider "acceptable." We could attempt to use an F-test of sorts, but with a multi-property weighted scoring function this gets more ambiguous.

Instead, I use the variation between the short, general, and long parameters as well as the scoring function to assign the following standard deviations:

eps_CH: 0.5 K
sig_CH: 0.05 A
eps_C: 0.1 K
sig_C: 0.08 A

The covariance between eps and sigma was assigned by visual inspection of the scoring function.

Here is my first pass at assigning uncertainties. Since the challenge compound is a "long" branched alkane, I am using the "long" parameter set as my maximum likelihood. Note that I use the same plot region as that of Mick et al. in an attempt to help visualize the scoring function compared with the MCMC points:

Does anyone object to these uncertainties? In other words, do you think they are too large, too small, or that the correlation between eps and sig is too strong or too weak?

If no one objects, I will begin simulating 200 different MCMC parameter sets which are independently sampled from the CH3, CH2, CH, and C parameter spaces.

Alternatively, if @jpotoff or @msoroush have the actual scoring function values, we could fit a model to this surface and perform MCMC on that model. This would certainly be more rigorous, but I was not sure if the scoring function values are readily accessible for all the different parameter sets.

@mostafa-razavi

We will probably want to use smaller system sizes for the larger compounds to reduce computational costs. We saw that for the ITIC approach we needed about 800 molecules for n-octane. However, most studies suggest system size effects are less pronounced for viscosity. Most studies use between 120-500 molecules. We will want to modify our script so that we can test for finite size effects. Hopefully we see similar results for 200 molecules as 800. Otherwise, our simulations will take much longer.

We are simulating the normal and branched alkanes to determine the adequacy of the united-atom Mie lambda-6 potential to predict viscosity. Although the ultimate goal is to predict high pressure viscosities, currently we are investigating the ability to predict saturation viscosity. In order to account for uncertainties in the parameter set (epsilon/sigma) for a given lambda, we have performed a Bayesian MCMC uncertainty analysis.

The figure below shows the feasible CH3 parameter sets for different values of lambda:

This figure shows the percent deviations in liquid density and vapor pressure for the different parameter sets:

As we can see, the 15-6 or 16-6 potential appear to be the best options for predicting VLE.

We then simulate a subset of these epsilon/sigma parameter sets to determine the overall uncertainty in viscosity for a given lambda.

Here are the viscosity results for ethane (just lambda =13-16), where the error bars represent the distribution in epsilon/sigma sets:

Clearly, the viscosity is over-predicted with the UA Mie potential, and the bias increases with increasing lambda. Although the 13-6 potential is better at predicting viscosity, recall that it is significantly worse at predicting VLE than the 15-6 or 16-6. In other words, there is a significant trade-off between accurate prediction of VLE and accurate prediction of viscosity. This suggests that we need an alternative (perhaps all-atom) model to accurately predict both VLE and viscosity.

We are repeating this analysis now for propane, butane, and n-octane where we include the uncertainty in the CH2 parameters.

@mrshirts

We cannot use MTTK because we are using fixed bonds (LINCS) . Also, we had been using velocity Verlet (md-vv) which is not compatible with Parrinello-Rahman. So we had two options in GROMACS, either run md-vv with the Berendsen barostat (which has been shown to not sample from the correct ensemble) or use the leapfrog integrator (md) with Parrinello-Rahman. For most cases, md is sufficiently accurate, so we probably want to use md with Parrinello-Rahman.

In order to support this decision, we have performed simulations using both barostats (note that currently the results for Parrinello-Rahman have tau-p of 5.0. We are rerunning these simulations with tau-p of 2 to see if that changes anything). The error bars from this study (red circles and blue squares) represent 95% confidence intervals. It is unclear what type of errors Nieto-Draghi reported, so we did not attempt to modify their uncertainties when plotting the error bars for the literature (green plus).

We found it interesting that the volume distributions have approximately the same width. Also, note that Parrinello-Rahman agrees more closely with the prescribed pressure for each state point except the lowest P (1 bar). Furthermore, the box volume and pressure show a systematic deviation from the literature values.

@mostafa-razavi

I had not asked you previously to run simulations for the challenge compound because it did not seem necessary and there was no clear distinction between which jobs you or I should run (i.e. we aren't simulating various force fields). However, I am nervous about having enough computational resources to finish all the required simulations by the challenge deadline (10/7). The main reason for this problem is that at the highest pressures we need extremely long simulations, and sometimes I have to rerun those simulations because they were not long enough. I cannot simply overshoot with the simulation time because the post data analysis is also very slow as the file size becomes inordinately large.

I am currently running simulations using 400 molecules, even though I know that a 200 molecule system would run about twice as fast and should not deviate considerably from the 400 molecule results. I would like to transition to a 200 molecule system, and I thought that it would make sense to include you in this effort.

Would you be able to simulate 224TMHexane with Potoff, 1.0 nm cut-off, 200 molecules, LINCS, 50-60 replicates, NPT first, for rho = 0 to 7?

I will run rho = 8 to 12 (the higher pressure, long simulations). I need to update the GitHub repository files, so don't start running any simulations until I give you the go ahead.

Please let me know if you have any concerns

Thanks

@mostafa-razavi @Kaiveria

Now that we have resolved the main issues in GreenKubo_analyze it is important for us to go back and reprocess the data. The three step process is simple:

1. Make sure you have the most recent version of GreenKubo_analyze.py
2. Edit the parameters in Rerun_Green_Kubo_Analyze.sh (e.g., Compound, Model, Conditions...) for the given system you want to reanalyze
3. Execute Rerun_Green_Kubo_Analyze.sh

The Rerun_Green_Kubo_Analyze.sh script currently has the parameters for rerunning IC8H18 for TAMie with N400. The error bars on this simulation have the problem that we have fixed in GreenKubo_analyze. So @mostafa-razavi could you reprocess this system and post your results this week? Eventually you will want to reprocess all of your results, but in the immediate future I just need this system for a presentation I am giving at the NIST TRC Consortium.

@mostafa-razavi

Currently my simulation set-up is the following:

Two energy minimizations (very fast)
Short equilibration (100-500 ps)
Full production period (1-4 ns)
Viscosity period (1 ns)

We probably do not need this long production period. In fact, we don't really need it at all if we are confident we have equilibrated well. I propose that we increase the equilibration period to about 1 ns for each simulation and remove the full production period altogether. This will have a large impact on the larger compounds, where we were going to be wasting a lot of CPU time on the production periods that are discarded.

@jrelliottoh

I validated that the pressure-viscosity coefficient is of a similar magnitude and follows a similar trend to those from the Bair models.

Here I reproduced the pressure plots provided by Bair:

Here I computed the pressure-viscosity coefficients:

@mostafa-razavi

@Kaiveria and I have validated our TraPPE results for n-octane by comparing the viscosity-pressure dependence with that reported in the literature from two different simulation studies. The study from Kioupis (and Ed Maginn) is from 2000 while Nieto-Draghi et al. (Philippe Ungerer and Bernard Rousseau) is from 2006.

Notably, Kioupis and Maginn used a harmonic bond potential, while Nieto-Draghi et al. appear to use fixed bonds (although they never say this explicitly, they did not report using a harmonic potential). For this reason, we performed two sets of simulations, one with fixed bonds and one with the same harmonic potential used by Kioupis and Maginn.

There are a few reasons why our results demonstrate a small discrepancy with the literature:

1. Nieto-Draghi reported their results at 347.52 K and purported that this was the same temperature as the Kioupis results. In fact, Kioupis and Maginn reported a value of 348.15 K. Furthermore, we accidentally performed our simulations at 347 K, which might explain why our pressures are slightly less than those reported by Nieto-Draghi, despite having the exact same density.

2. Also, note that Kioupis and Maginn used nonequilibrium molecular dynamics (NEMD) while Nieto-Draghi et al. used EMD with the Einstein relationship at short times. The extrapolation of NEMD results to the zero shear limit using the Carreau model typically over predicts viscosity. The Einstein relation at short times is a cheap, but not as rigorous approach, for obtaining reliable viscosities.

For these reasons, we have deemed our results to be sufficiently consistent with those from the literature. However, we will rerun the simulations at 347.52 K to see if this resolves the small shift in pressure. Furthermore, we have some simulations running that start with the NPT ensemble to determine the average density. These new simulations use the Kioupis pressure to initialize the system whereas the old simulations (shown above) used the Nieto-Draghi density as the initial configuration.

These results will be included as supporting information for validation.

This is a list of the simulations we need to run

Put a * next to high priority tasks
Put your initials next to a task that you are currently working on

Saturation viscosity

New list of compounds

C12H26

• Potoff (SMR: done)

• TraPPE

• TAMie

C16H34 (not working)

• Potoff

• TraPPE

• TAMie

C22H46

• Potoff

• TraPPE (SMR: running)

• TAMie

22DM-butane (probably won't run this)

• Potoff

• TraPPE

23DM-butane

• Potoff

• TraPPE

• TAMie

3-methylpentane

• Potoff

• TraPPE

• TAMie

Old list of compounds

C2H6

• Potoff

• TraPPE

• TAMie

• Exp-6

• TraPPE-2

• MCMC 12-6

• MCMC 13-6

• MCMC 14-6

• MCMC 15-6

• MCMC 16-6

C3H8

• Potoff

• TraPPE

• TAMie

• Exp-6

• MCMC 14-6 (low priority)

• MCMC 16-6

C4H10

• Potoff

• TraPPE

• TAMie

• Exp-6

• MCMC 14-6 (low priority)

• MCMC 16-6

C8H18

• Potoff

• TraPPE

• TAMie

• Exp-6

• MCMC 14-6 (low priority)

• MCMC 16-6

IC4H10

• Potoff

• TraPPE

• TAMie (SMR: Results are very similar to TraPPE)

• AUA4

IC5H12

• Potoff

• TraPPE

• TAMie

IC6H14

• Potoff (SMR: Done. Ran again due to wrong previous Potoff parameters for IC6)

• TraPPE

• TAMie

IC8H18

• Potoff

• TraPPE

NEOC5H12

• Potoff

• TraPPE

• AUA4

Finite size effects

C3H8, Potoff:

• N100
• N200
• N400
• N800

C4H10, Potoff:

• N100
• N200
• N400
• N800

C4H10, TraPPE:

• N200
• N400
• N800

C8H18, TraPPE:

• N200
• N400
• N800

IC8H18, TraPPE:

• N200
• N400
• N800

Simulation length

C8H18, TraPPE, N400

Saturation

• 0.5 ns
• 1 ns
• 2 ns
• 4 ns

High pressure

• 1 ns
• 2 ns
• 4 ns

Harmonic bonds

TraPPE, N400

• C3H8

• C4H10

• C8H18

High pressure

293 K

TraPPE, N400

• C3H8

• C4H10

• C8H18

• IC5H12

• 3MPentane (SMR: Running)

Potoff, N400

• C3H8

• C4H10

• C8H18

TAMie, N400 (high priority)

• C3H8

• C4H10

• C8H18

MCMC Mie 16-6, N400

• C3H8

• C4H10

• C8H18

MCMC Mie 14-6, N400

• C3H8

• C4H10

• C8H18

347 K

TraPPE, N400, C8H18 (Validation)

• Harmonic

• LINCS

@Kaiveria @mostafa-razavi @jrelliottoh

OK, now that we have finalized our manuscript for the special issue of FPE, it is time to move our attention to the challenge compound. The current plan is to simulate 224TMH using the Potoff model without modifying any of the parameters. We can then consider refining the Potoff parameters using MBAR-GCMC or MBAR-ITIC. Also, I have begun transferring the Potoff model into an ex-LJ potential (e.g., 11-10-9-8-7-6) which might perform better at extreme pressures.

@jpotoff @msoroush @mostafa-razavi @jrelliottoh

The viscosity at high pressures is quite sensitive to the torsional model. I know that we discussed how we do not want to modify the torsional parameters to match viscosity data because we still want to match ab initio values. I should mention that there is some debate as to whether these ab initio values (computed with a low level of theory and small basis set) should still be used for united-atom models. That being said, I agree that we do not want to reoptimize the torsional parameters. However, I do think we should assign some degree of uncertainty to the torsional potential.

Unless we already have some ab initio data we want to use to rigorously quantify the uncertainty, we can just assign some reasonable percent uncertainty. For example, if we believe our torsional potential is reliable to within 10%, we could assign this type of uncertainty:

We can use MCMC to randomly sample the parameter sets that meet this level of uncertainty.

Whereas the uncertainty in the Mie 16-6 parameters is somewhat negligible compared to the inherent simulation uncertainty in viscosity, the torsional parameter uncertainty might play a larger role.

@jrelliottoh @mostafa-razavi @Kaiveria

As the highest pressure state points require much longer simulations (48 ns), we are running into memory limits for the post-processing. Essentially, we can only process a few simulations at a time due to the file size, and these files take a day to process. This approach will take too long with the deadline in a few weeks. I have tried running simulations with less frequent outputs, but I have seen mixed results. In some cases, the reduced frequency seemed to have little impact, while in other cases the plateau was shifted to a higher viscosity. My initial attempt was reducing the output frequency by a factor of 6. Maybe that was just too dramatic. I will now try reducing it by just a factor of 2, 3, or 4.

Email from @mostafa-razavi :

When I run GreenKubo_analyze.py I get this error:

mostafa@mostafa:~/Git/IFPSC_10/C4H10/Gromacs/Mie_16_Bayesian_Viscosity_200_MCMC$ python ~/Git/IFPSC_10/Scripts/GreenKubo_analyze.py

Traceback (most recent call last):
File "/home/mostafa/Git/IFPSC_10/Scripts/GreenKubo_analyze.py", line 1343, in
main()
File "/home/mostafa/Git/IFPSC_10/Scripts/GreenKubo_analyze.py", line 1330, in main
MCMC_Mie = GreenKubo_MCMC(args.ilow,args.ihigh)
File "/home/mostafa/Git/IFPSC_10/Scripts/GreenKubo_analyze.py", line 53, in init
self.base_path_all = self.gen_base_path()
File "/home/mostafa/Git/IFPSC_10/Scripts/GreenKubo_analyze.py", line 379, in gen_base_path
for iMCMC in np.arange(ilow,ihigh+1):
TypeError: unsupported operand type(s) for +: 'NoneType' and 'int'

I have a few questions:

1- Does the current version of GreenKubo_analyze.py work for n-butane as it is? If not, what changes should I make?

2- How (using what command-line arguments) should I run GreenKubo_analyze.py? What directory should I run this script in?

@mostafa-razavi
I just created a Mie_14_Bayesian_Viscosity_200_MCMC directory for butane. I modified the nAlkanes... script so that model is Mie_14... and compound is C4H10. Can you start those simulations? It should take the same length of time as your last simulations, i.e. about 6 days.

I am currently rerunning the propane Mie_16 200 MCMC results (it takes about 3 days). The first time I accidentally used the box sizes from ethane, so they were highly compressed liquid simulations. On the bright side, I was able to get meaningful viscosities of around 20 cP using Green-Kubo. Most people would argue you need NEMD for that high of a viscosity, but with replicates you can use EMD.

The Potoff viscosities for the challenge compound at 1000 MPa are much greater than we anticipated. Even a 24 ns run is not long enough to observe a clear plateau:

We can try running for even longer, but at some point I am worried about the disc space. A 24 ns run takes up about 0.82 GB, while a 16 ns run takes up about 0.54 GB. So a 30 ns run should take up about 1 GB per replicate. But I think we need to try this.

@Kaiveria

Do you have any concerns with this? Other than how slow the post simulation analysis will be.

Here are some scripts that I will post this week:

nAlkane script for the Lennard-Jones native potential
nAlkane script for literature force fields

branchedAlkane scripts for tabulated, Lennard-Jones, MCMC, literature, etc.

We will generate radial distribution functions to help explain why the Potoff Mie potential behaves worse at higher densities. Specifically, we are interested in how the initial rise of the RDF changes as a function of density, e.g. the distance where RDF > 0, the slope of RDF with respect to r, the distance where RDF = RDF_max.

@jrelliottoh @mostafa-razavi @Kaiveria

It became evident that it was not an ideal plan to have the manuscript in the same repository as all the scripts, etc. So the manuscript has been moved to this GitHub repo:

https://github.com/ramess101/Special_Issue_FPE