United-atom Mie lambda-6 viscosities

[ramess101]

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.

[mostafa-razavi]

@ramess101 Results of Mie n-butane when n=16 is attached. Do we have a script to generate the above plot, i.e. the viscosity deviation from Refprop vs temperature? I also attach the text files containing the data:

Results_Mie_16_Bayesian_Viscosity_200_MCMC.zip

[ramess101]

@mostafa-razavi

Yes, I have a script but I will clean it up before I post it.

[ramess101]

@mostafa-razavi

OK, I have uploaded "compare_TDE_REFPROP.py" that compares the simulations with the REFPROP correlation and TDE experimental data. The header for this file explains what files you need to call it. I have also made it so that it will be called automatically in the future after GreenKubo_analyze.py.

Note that I have created new files called "nAlkanes_SaturatedViscosity" and "branchedAlkanes_SaturatedViscosity". These files will use the native Lennard-Jones function in Gromacs if the force field is a 12-6 potential (this is about twice as fast as tabulated potential). Also, note that I also renamed some directories. Specifically, "$Compound"/Gromacs/Tabulated is now just "$Compound"/Gromacs/Gromacs_input because the .gro and .top files are the same for tabulated and Lennard-Jones.

[ramess101]

@mostafa-razavi

This is what you should get with the files you sent me for 16-6 MCMC:

Note that the large uncertainty at low Tsat appears to be a bug in my bootstrap analysis code. We will need to work on that, but it does not change the primary conclusion that the Mie 16-6 potential actually under predicts viscosity for butane!

This is consistent with what I get for 16-6 MCMC for propane:

And for propane with Potoff:

Note that the abnormal value at 210 K is likely because of data analysis as well. We will figure that out, but for now it is good enough.

This is very curious that the 16-6 is actually under predicting viscosity.

The 14-6 (TAMie, which is also an AUA model) appears to under predict viscosity to an even higher degree:

This explains, in part, why Gordon reported a 20-6 CH2 potential for n-alkanes when trying to reproduce viscosity data.

[ramess101]

@mostafa-razavi

OK, so this under prediction is actually fairly well-documented. It is a result of the united-atom approach. Gordon argued that you can change the repulsive exponent to overcome this issue, but his VLE results were worse. We have shown that for adequate VLE the viscosity will still be under predicted. Ethane appears to be an exception.

Gordon provided some results for n-alkanes using the TraPPE model (open circles):

I am currently running TraPPE for propane, but it would be good to get TraPPE for butane finished so we can make a comparison. I will start TraPPE for butane unless your Potoff runs finish before mine, since you have done all the other butane runs.

One important question is the impact of fixed bonds. We are using constraints while Gordon used LAMMPS, so his larger compounds could not be constrained and he never explicitly says if his smaller compounds are constrained or what bond potential he is using.

[mostafa-razavi]

@ramess101 Here is the Potoff butane three lower densities that finished. I'll run the other two densities again to complete this.

@ramess101 I updated the plot to include rho3 and rho4.

Note that the large uncertainty at low Tsat appears to be a bug in my bootstrap analysis code. We will need to work on that, but it does not change the primary conclusion that the Mie 16-6 potential actually under predicts viscosity for butane!

Is this problem solved?

[ramess101]

@mostafa-razavi

Those results look good.

Is this problem solved?

No. Sorry, not high on my priority list. But it should be.

[ramess101]

@mostafa-razavi

Try using "GreenKubo_analyze" now to analyze the low density Potoff butane run. I merged the GreenKubo_analyze_fitw8.py code and modified it to run faster. This is only a minor fix, but if this doesn't improve the results I have a couple other things I can try.

[ramess101]

@mostafa-razavi

I think that it might be important to look at the viscosity trend approaching the triple point. For example, while the Potoff model does fairly well at these moderate viscosities (between 0.1 and 1 cP), I am curious what it looks like for 1 to 10 cP. I was thinking about including simulation conditions for propane at T= 90, 100, 110, 120, and 130 K. Recall that these were the Potoff propane results between 160 and 320 K:

What do you think?

[ramess101]

@mostafa-razavi

Here are the Potoff propane results at lower T. I am surprised that it works as well as it does. I should probably use more replicates since the uncertainties get quite large.

However, recall that these simulations are performed at the REFPROP Tsat and rholsat. I am assuming that each force field gives a reasonable rholsat.

[ramess101]

@mostafa-razavi

The N=200 octane TraPPE results are finished.

It will be interesting to see how much better the Potoff model does. Note that I only used 30 replicates at the lowest temperature, while the rest had 40. This can partially explain the large uncertainty. Also, using a larger system size should decrease the uncertainties. But the real reason for these large uncertainties is still the GreenKubo analysis. It is really not a bug, we just need to discuss how we want to compute the uncertainties. Right now, it allows for the cut-off time to be a random variable, which can occasionally lead to abnormal viscosity values.

For example, this is the distribution of bootstrapped viscosities at the two lowest T:

T=256 K:

T = 321 K

Clearly, most of the values are in the higher viscosity region for T=256 K. But, there is a second mode at lower viscosities. This mode is very pronounced for T=321 K. If we used a constant cut-off this would likely not be the case (I can test this quite easily).

For a comparison, this is what the distributions look like at the higher temperatures:

[ramess101]

@mostafa-razavi

I created a script called compare_force_fields that produces a plot for multiple compounds and force fields. These are the results for n-alkanes that I have compiled so far. It is incomplete, but already starting to paint a clear picture that Potoff does a lot better than TraPPE and slightly better than TAMie:

[ramess101]

@mostafa-razavi

I performed some non-saturation simulations with the Potoff potential for propane. As expected, the viscosity is over predicted for the higher densities. However, if we plot the viscosity with respect to pressure it actually agrees very nicely. Note that the REFPROP correlation is only valid to 1000 MPa, so we do not really know if that highest point is good or not. I will run some simulations below 1000 MPa to see if the Potoff model is actually reliable enough for this case.

[ramess101]

@mostafa-razavi

I performed four more points with Potoff propane at high densities:

Clearly something is wrong with the uncertainty for the 1000 MPa point, but it is rather remarkable that the Mie 16-6 potential works so well up to 700 MPa.

I now have the TraPPE propane runs going.

[ramess101]

@mostafa-razavi

These are the TraPPE propane results:

So the Potoff model is much better at predicting even high pressure viscosities. However, both models are deficient when plotted with respect to density.

[ramess101]

@mostafa-razavi

I have a new GreenKubo_analyze code that is more robust in assigning uncertainties. For example, here I have rerun the bootstrap uncertainty at P around 700 and 1000 MPa, which previously had very large uncertainties:

I am currently reprocessing the other points. Moving forward make sure to use the new GreenKubo_analyze script.

[mostafa-razavi]

@ramess101

Potoff-C4-N200

Potoff-C4-N200.zip

[mostafa-razavi]

@ramess101

Potoff-C4-N800

Potoff-C4-N800.zip

[mostafa-razavi]

@ramess101

TraPPE-C4-N400

TraPPE-C4-N400.zip

[mostafa-razavi]

@ramess101

TAMie-C4-N400

TAMie-C4-N400.zip

[mostafa-razavi]

@ramess101

TraPPE-IC4H10-N800

TraPPE-IC4H10-N800.zip

[ramess101]

@mostafa-razavi

Wow, those uncertainties are very small for TraPPE IC4H10 N800 (compared to what we usually see where the uncertainy is +- 5 %). I would assume this is because of the larger system size, but previously in Issue #5 the N800 uncertainties were not consistently smaller. The results in Issue #5 used the old GreenKubo_analyze, so maybe now with the better GreenKubo_analyze the uncertainties are consistently smaller for larger systems.

**Correction, since viscosity is a collective property the uncertainties do not depend on the number of molecules.

[mostafa-razavi]

TraPPE-IC5H12-N400

@ramess101 There seems to be a problem in TDE data, namely there are two trends that are significantly different. What do you think?

TraPPE-IC5H12-N400.zip

[mostafa-razavi]

TraPPE-IC6H14-N400

@ramess101 TraPPE-IC6H14-N400.zip

[ramess101]

@mostafa-razavi

There seems to be a problem in TDE data, namely there are two trends that are significantly different. What do you think?

I trust the TDE data that agree closely with the REFPROP correlation. The % error for TraPPE appears to be consistent, namely, 20 to 50 % under prediction that increases at lower temperatures.

[ramess101]

@mostafa-razavi

Here are the n-butane high pressure results. Note that the REFPROP correlation is only recommended up to 200 MPa (700 kg/m3).

Potoff:

TraPPE:

Again, we see that the Potoff potential over predicts viscosity at high densities. Because of the large uncertainties in REFPROP above 200 MPa, it is difficult to ascertain whether or not the Potoff model is unreliable at high pressures. Clearly, the TraPPE model under predict viscosities at all conditions, and appears to get slightly worse at higher pressures.

[mostafa-razavi]

TAMie-C8H18-N800

@ramess101 TAMie-C8-N800.zip

(Updated using Rerun_Green_Kubo script) TAMie-C8H18-N800.zip

[mostafa-razavi]

TraPPE-IC8H18-N400

@ramess101 TraPPE-IC8H18-N400.zip

[mostafa-razavi]

TAMie-IC4H10-N400?

@ramess101 Could you confirm if ran the right script for this simulation? I would expect TAMie to do better than this. Right?

TAMie-IC4H10-N400.zip

Update: (2/7/2018)

TAMie-IC4H10-N400

@ramess101 @Kaiveria Here is what I get from the new simulation. It's very close to TraPPE but slightly different. You might want to run one of these points on your computer to see if my simulation is correct or not. compare_TDE_REFPROP_sat.pdf TAMie-IC4H10-N400.zip

[ramess101]

@mostafa-razavi

Those results look almost identical to the TraPPE results. I can take a closer look, but I am pretty confident these are not the TAMie results.

[ramess101]

@mostafa-razavi

Your "IC4H10_job_2018_05_23_18_37_17" file says that you submitted TAMie. However, the only way to be certain is to find the IC4H10.top file, which should be at IC4H10/Gromacs/Saturation_Viscosity/TAMie_N400/MCMC_x/IC4H10.top. Can you post that file?

Also, these are the parameters you should have:

[ atomtypes ] ; _O indicates the (pseduo)atom is bonded to an oxygen. _C indicates the ; (pseudo)atom is bonded to a carbon. ; name bond_type mass charge ptype C6 Clam CH3_C CH3 15.03450 0.000 A 0.00819752345327 2.33014239486e-06; [CH3]-CHX CH2_C CH2 14.02660 0.000 A 0.00631986287433 4.48496002146e-06; (CHX)2-[CH2] CH_C CH 13.01900 0.000 A 0.002764768497 3.64762596488e-06; (CHX)3-[CH] C_C C 12.01100 0.000 A 0.002764768497 3.64762596488e-06; (CHX)4-[C]

[ nonbond_params ] ; i j func C6 Clam CH3_C CH2_C 1 0.00726864063375 3.3075463302e-06; cross-interactions for CH3_C and CH2_C CH3_C CH_C 1 0.00489378022934 3.10910301117e-06; cross-interactions for CH3_C and CH_C CH3_C C_C 1 0.00489378022934 3.10910301117e-06; cross-interactions for CH3_C and C_C CH2_C CH_C 1 0.00419893969751 4.08742151221e-06; cross-interactions for CH2_C and CH_C CH2_C C_C 1 0.00419893969751 4.08742151221e-06; cross-interactions for CH2_C and C_C CH_C C_C 1 0.00419893969751 4.08742151221e-06; cross-interactions for CH_C and C_C

[mostafa-razavi]

@ramess101 IC4H10.top.txt is one of the files you asked for. It seems to be identical to the values you posted above. Could there be a problem with the compare_TDE_REFPROP.py.txt script that I used? I used this old version because, the original compare_TDE_REFPROP.py.txt script in the "Script" folder give me an error about not having a value for MW for that compound in

Mw_list = {'C3H8':44.096,'C4H10':58.122}

[ramess101]

@mostafa-razavi

I don't think the compare_TDE_REFPROP.py file is the problem, because the eta_fit_rho0, etc. figures also look like the TraPPE results. Is it possible that you called GreenKubo_analyze from the wrong directory somehow? The script should do this automatically, but it could be an issue if you had to call it manually for some reason.

[mostafa-razavi]

TraPPE-IC8H18-N200

TraPPE-IC8H18-N200.zip

[mostafa-razavi]

TAMie-IC5H12-N400

TAMie-IC5H12-N400.zip

[ramess101]

This plot compares Potoff, TraPPE, and TAMie for propane at 293 K and high pressures:

[mostafa-razavi]

TAMie-IC6H14-N400

TAMie-IC6H14-N400.zip

[ramess101]

@mostafa-razavi

These all look really good. I am still confused about TAMie-IC4H10 that looked almost exactly like TraPPE-IC4H10. We need to figure that out at some point, or we can just rerun it sometime.

Mostafa, for your next simulations, could you run C4H10 and C8H18 at high pressures (T293highP) for the TAMie force field? You can use the "nAlkanes_T293highPViscosity" file. (We have some better scripts for branched alkanes that are more robust, but we are still testing to make sure they work well.) Right now the nAlkanes_T293highPViscosity script is set up to run C4H10 with TAMie for the 5 T293highP state points, so you shouldn't need to change much. Let me know if you have any issues with the script though.

[ramess101]

@mostafa-razavi

I just wanted to see if you were able to get the T293highP viscosities running for C4H10 and C8H18 with TAMie. I am currently out of town, so please ask @Kaiveria if you have any questions with the scripts. Thanks

[mostafa-razavi]

TraPPE-IC8H18-N800

TraPPE-IC8H18-N800.zip

[mostafa-razavi]

I just wanted to see if you were able to get the T293highP viscosities running for C4H10 and C8H18 with TAMie. I am currently out of town, so please ask @Kaiveria if you have any questions with the scripts. Thanks

@ramess101 Right now I'm running TAMie_C4H10_N400_T293highP. One of the computers that I'm using is under maintenance and I do not have access to it. I will check if I can use our other computers to run TAMie_C8H18_N400_T293highP.

[mostafa-razavi]

TAMie-C4H10-N400-T293highP

compare_REFPROP_T293highP.pdf TAMie_C4H10_N400_T293highP.zip

[mostafa-razavi]

TAMie_C8H18_N400_T293highP

@ramess101 Do we have any experimental data to complete the following plot? compare_REFPROP_T293highP.pdf TAMie_C8H18_N400_T293highP.zip

[ramess101]

@mostafa-razavi

Do we have any experimental data to complete the following plot?

Yes we do. When I create the final figures the data will be included.

I think we have completed all of the simulations that are currently on the work plan Issue #9 (at least those that really need to be run). We have a few more molecules that should be added, so check back on Issue #9 by tomorrow for the next batch to run. Also, it is very important that you double check that you have the most recent .top files. We just discovered a typo that impacts the larger branched alkanes (IC5, IC6, IC8).

[ramess101]

@mostafa-razavi

I just realized that you are already running the Mie 16-6 MCMC for n-octane. This is probably not the best use of computational resources. I don't think we are actually going to include the Mie MCMC results in our manuscript. So I would probably kill them (unless they are nearly finished) once we have selected the next system to run.

[mostafa-razavi]

TraPPE-C12H26-N400

(Number of Replicates: 48, Duration: 2 days and 6 hrs, on 24 cores)

compare_TDE_REFPROP_sat.pdf TraPPE-C12H26-N400.zip

[mostafa-razavi]

TraPPE-3MPentane-N400

compare_TDE_REFPROP_sat.pdf TraPPE-3MPentane-N400.zip

[mostafa-razavi]

TraPPE-IC5H12-N400

compare_REFPROP_T293highP.pdf TraPPE-IC5H12-N400.zip

[mostafa-razavi]

TraPPE-3MPentane-N400-T293highP

compare_REFPROP_T293highP.pdf TraPPE-3MPentane-N400-T293highP.ZIP

[mostafa-razavi]

Potoff-IC6H14-N400

compare_TDE_REFPROP_sat.pdf Potoff-IC6H14-N400.zip