Special Issue of FPE

[ramess101]

From Mostafa:

Richard,

Thank you for forwarding this email. This special issue sounds like a good outlet for our research. Any ideas what could the title of our submission be?

Generally, I can imagine a few routes for our research:

1- A comparison between different potential models: shifted-force, force-switch, extended LJ, Mie, Yukawa, etc. 2- Developing force fields for additional compounds using MBAR/ITIC 3- Building an online framework for customized force field development similar to SpeadMD where people can develop their own force fields by taking advantage of MBAR 4- Developing a transport-properties optimized force field: participating in IFPSC challenge (See http://fluidproperties.org/10th) could further motivate us for focusing on transport properties

Obviously, these are brainstorming ideas and their feasibility and priority can be discussed. It would be nice if we complete such a list and prioritize them.

[ramess101]

Mostafa,

These are some good ideas. Here are my thoughts regarding each point.

1. A comparison between different potentials would certainly be an interesting study. We would need to determine what exactly we want to compare. A comparison of how well they predict VLE properties, for example, really depends on an adequate parameterization of the potential, which depends on the objective function and a good optimizer. Was that what you were envisioning? With MBAR/ITIC this should be possible, but I think it would be a lot of work to compare several functional forms that have all been parameterized in the same way for a wide range of compounds.
2. This is our goal for collaborating with Jeffrey Potoff. However, I don’t think we will have anything ready by August.
3. I am somewhat hesitant to go down this path. It seems like a lot of work and I am not confident in my abilities to build an online framework. Also, I think we probably want our method to gain some traction before we have any hope for people actually seeking out the tool to utilize it. For now, I think having useable script files on GitHub is the best option.
4. I actually am already working on. I have a summer intern that will be working with me from May-August on an IFPSC 10 submission. I would welcome your contributions, if you would like to join our efforts. I have also discussed this briefly with two other groups (Ed Maginn and DIPPR). Currently I am working on developing a force field that has the correct PVT behavior at high pressures. I have an upcoming publication which demonstrates that the Mie potential is not adequate for predicting PVT at high pressures as it systematically under-predicts the density (i.e. over-predicts the pressure). Instead, we need to either use an anisotropic-united-atom or an all-atom model. I am working on developing an AUA LJ 12-6 model because the TraPPE-2 model showed promising results. As for the publication, I believe Fluid Phase Equilibria has a separate edition specifically for IFPSC submissions, right? So I don’t think we would submit our results to the August FPE edition. At least that was my understanding.

I was thinking that we could do a combination of your ideas 1 and 2. Specifically, we could demonstrate the speed of MBAR/ITIC to reoptimize parameters for a force-switch potential. For example, we could publish our results where we reoptimized the TraPPE-2 force field using a force-switch. We could then repeat this for a large set of pseudo-diatomics, i.e. two UA site molecules. This would be a good starting point for developing a larger set of force field parameters. However, I do not want to commit too much of my time to this endeavor since I envision IFPSC taking a lot of my time already. Perhaps we could just publish the TraPPE-2 force-switch results that we already have. This would demonstrate “advances in algorithms for parameter optimization” and would raise the issue of force-switch potentials. The title would be something like: Rapid conversion of the non-bonded parameters from a cut-off to a force-switch potential. What do you think?

Rich

[ramess101]

Mostafa,

I wanted to follow-up on this conversation. I was thinking about the work we did converting the TraPPE-2 to a force switch for ethane. Our methodology was to reoptimize epsilon and sigma based on rhol and Pvsat. The “optimal” parameter set resulted in accurate rhol and Pvsat, but not identical to TraPPE-2, i.e. we were slightly better in rhol but slightly worse in Pvsat. This is because our objective function was inherently different than that used by Siepmann, especially since ITIC is testing a different region of temperatures. This was one of my primary concerns for optimizing other compounds where TraPPE does not accurately predict Pvsat. However, if our goal is to instead demonstrate how MBAR can convert an existing cut-off force field to a force-switch, I think we could use a completely different approach that should be a lot simpler and would allow us to really call this the TraPPE-fs or Potoff-fs, etc.

Rather than using ITIC to obtain rhol and Pvsat and then defining our own objective function, we can determine the parameter set that reproduces the energy and pressure of the original force field. In other words, we simulate a few state points (could be all of the ITIC conditions but not necessarily) and compute the energy and pressure for the original cut-off based force field, with the corresponding uncertainties. We then perform a weighted optimization of U and P, where the weights are the inverse variance in U and P, to find the optimal epsilon and sigma for the force-switch force field.

This approach should be a lot easier/faster to implement since it would not matter if the original force field is accurate for rhol or Pvsat. In fact, we could start on this right away, since we already have the TraPPE simulations for n-alkanes and branched alkanes.

I think a study like this would be great for the FPE special issue. Also, this would be a logical preliminary step to optimizing the Mie-fs potential with Potoff. Basically, we can optimize a cut-off potential and then rapidly convert it to a force-switch potential. This methodology would be incredibly useful in practice for researchers that want to use a force-switch potential and do not know how to convert their cut-off based force field of interest.

Let me know what you think

Rich

[ramess101]

Richard,

Here is my response to points 1-4 and your follow-up email:

1- I agree that what you described requires a lot of effort, and we might not have much to present in August. I was also envisioning the fundamental questions such as why do we insist on force-switch and shifted-force potentials, and why should people use these potentials and not readily available TraPPE-UA. Here we could, for example, compare different truncation strategies in terms of conservation of energy using long NVE runs, and hopefully show that force-smoothed potentials have an advantage. Other questions that we could explore are similar to what you described in your point 4: "...the Mie potential is not adequate for predicting PVT at high pressures...". For example, which three-parameter potential model (such as Mie, extended LJ, etc.) gives us more flexibility for calculating Psat, rhoL, compressed liquid density, viscosity, etc.

2- When do you think we are ready to start that endeavor? How many papers are we expecting to be published from our collaboration with Jeffery Potoff? What would the tentative titles be?

3- Since I'm also doing a MS thesis in computer science, I appreciate the kind of exposure that I will have as a result of working on such projects. I do not think building an online framework is what we should be working on right now, but I do believe that after we gained some confidence about the models that we are developing, such framework will promote our work. We can decide when we reached to that point.

4- I am indeed interested in contributing to your IFPSC 10 submission. How are you planning to approach the problem? Are you thinking about developing an AUA or all-atom model for 2,2,4 Trimethylhexane?

5 (follow-up)- I have thought about this approach before and I briefly mentioned something similar in my MS thesis. I was able to quickly obtain a surface response in terms of epsilon and sigma of ethane using P and U of a few ITIC state points. I believe this might work if we carefully choose the best NVT points. Using MBAR, we can easily consider all ITIC points plus the points required for B2 calculation. This approach is especially useful for compounds for which we do not have experimental P and U data.

Since the time is somehow limited, I believe we can split tasks, i.e. I will concentrate more on FPE special issue, while you focus more on IFPSC. This way we will have two projects progressing at the same time. What do you think?

S. Mostafa Razavi

[ramess101]

Mostafa,

1 - OK, something like that sounds a lot more feasible. Do you already have some preliminary results to show the need for force-switch and shifted-force? I think it could be good to demonstrate or at least discuss the situations when a cut-off could have problems. Or, like you mentioned, we could use this as an opportunity to look at the extended LJ. Specifically, we could see if extended LJ extrapolates better to high pressures or to transport properties when parameterized with VLE data. Alternatively, we could test the Mie potentials in the literature (Potoff and TAMie) and see how well they do for viscosity. This would be a pretty simple study but I think it would meet the expectations of FPE.

2 - I think we won't have any publications for quite a while, maybe 9 months. We have several infrastructure issues to work out before we can go into production mode. And things are moving slow right now since this is more of a secondary activity.

3- Oh, I did not know you were getting an MS in computer science. I agree that this would give our work good exposure. So when we get our feet under us we can look into this.

4- Great! After seeing that the Mie potential will have serious issues at high pressures, we are planning on using an AUA model. Basically, we are going to use the TraPPE-2 interaction displacement distance for the CH3 sites. We will have all other sites at the carbons for simplicity and because the CH2 and CH sites are pretty well shielded. We are going to use MBAR to optimize the force field parameters by matching U and P of REFPROP to high pressures for isooctane, which is the best surrogate for 2,2,4-trimethylhexane. If the AUA model does not accurately predict viscosity for isooctane, we might consider using an all-atom model. But we would like to avoid that if possible.

5- Let me try to clarify exactly what I have proposed to see if indeed we are talking about the same thing.

These would be the steps to convert TraPPE-UA to TraPPE-UAfs:

1. Simulate the TraPPE-UA force field at some selected state points
2. Calculate the TraPPE-UA U/P at each state point
3. Optimize the epsilon and sigma for TraPPE-UAfs such that we minimize the deviation between the TraPPE-UA and TraPPE-UAfs U/P at each state point, NOT with respect to experimental (i.e. REFPROP) values

So we would not need to compute B2, or even rhol/Psat, and we do not even need experimental U/P data, since the goal is not really to improve our prediction of VLE or any other property. It is simply to provide a rapid conversion from a cut-off to a force-switch potential that is consistent with the original parameters.

6 - (FPE/IFPSC) OK, I think it is a good plan for me to focus on IFPSC and you can start working on the FPE special issue. But I think we are still not 100% sure what we want to do for FPE. We discussed three different options in Point 1 above and one more in Point 5. Specifically I think these are the options we are considering:

Option a) Discussing when force-switch or shifted-force are required/better than cut-off

Option b) Determine if the extended LJ is better for predicting VLE and high pressures (or viscosities) than the Mie potential

Option c) Report viscosity trends for Mie potentials

Option d) Demonstrate how MBAR can rapidly convert a cut-off potential to force-switch (or shifted-force, etc.)

I can see each of these being beneficial. I think the key questions to determine which option to choose are:

Which one would be the most important/impactful?

I think Options b and c. Options b and c are very practical and are data-driven. Although Options a and d are of great importance to us, they might not be important to the majority of practitioners.

Which one would lead to future publications, especially with the Potoff collaboration?

I think Options b and d. If the results for Option b are promising this could open a whole new line of force fields. Option d could be a good preliminary publication for the Mie-fs force field with Potoff.

For which one do we already have some preliminary results?

I know I have some for Options c and d. Do you have any results for Options a and b?

For which one are we the most confident that we can have something by August?

I think Options a, c, and d. I am afraid Option b could be more challenging than we anticipate.

Let me know what you think. Also let me know if you have any other options and questions to consider.

Rich

[ramess101]

Rich,

1- I do not have any preliminary results to show the need for force-switch and shifted-force.

2,3- OK.

4- Sounds like a good plan. Are we going to use a force-smoothed potential or tail correction? Do we have Refprop data for 2,2,4-trimethylhexane or 2,2,4-trimethylpentane? NIST WebBook does not have fluid properties data for these two compounds.

5- I think we are on the same page. I meant to say even though we do not directly use Psat in our optimization it is useful to compare the supercritical low density NVT points (i.e. the points required to calculate B2 value at suprcritical isotherm) between TraPPE-UA and TraPPE-UAfs and try to minimize the deviation between the two, because I believe those points have a great impact on vapor pressure and we want TraPPE-UA and TraPPE-UAfs to predict the same Psat behavior. Do you agree? Also, by mentioning Refprop I meant this whole idea of translating (converting cutoff to force-switch or shifted-force) force fields has one particular advantage: we can obtain a shifted-force or force-switch force field for a compound for which we do not have any Psat/rhoL experimental data by translating an existing force field (e.g. TraPPE-UA).

6- I agree with your answers to the the raised questions. Since I do not have data for options a and b (I do have ITIC results for 15 TraPPE-UA compounds which is a little in favor of option d), it seems like options c and d are winning. Option c is more in line with IFPSC and I think our efforts will complement each other and the chances of meeting the deadlines will improve. On the other hand, as you said, option d seems like a good step towards our long-term goals and Potoff collaboration. I will discuss this with Dr. Elliott today and I'll let you know.

S. Mostafa Razavi

[ramess101]

Mostafa,

1- OK. We should at least have some idea as to when we would need/want a force-switch over a cut-off.

4- I was planning on using cut-off with tail corrections because I know that the basis functions approach works. We can always convert this to a force switch afterwards. 2,2,4-trimethylhexane is the compound of interest, and relatively little data are available for it. However, 2,2,4-trimethylpentane, a.k.a. isooctane, is a very well studied compound. Isooctane is included in REFPROP and we even have viscosity data to fairly high pressures. So that is why we are going to test our force field on 2,2,4-trimethylpentane before using it on 2,2,4-trimethylhexane.

5- OK, that makes sense. Yes, I think we are on the same page. The low density values are certainly important for accurate Pvsat. Since they typically have larger uncertainties they might not have a large weighting in our objective function. We will need to figure out how important they really are. I would think that if we match Z-1/rho on IT at high densities, Z on IC, and U on IC we probably will match Z-1/rho at low densities anyways since we are not changing the potential significantly, i.e. it is still a LJ 12-6. But we would want to test this. Yes, one of the key aspects of this project is that we do not need any experimental data, since we are simply matching an existing force field.

6- I agree that Options c) and d) probably make the most sense.

6c- Option c) would line up very nicely with IFPSC. I think we could test the viscosity at high pressures using the Potoff force field for the normal and branched alkanes in REFPROP. I have two concerns with Option c) though. First, I do not have the computational resources to perform a large amount of simulations, which would be required to have precise estimates of viscosity for a large number of compounds. Second, we would need to use GROMACS to have fixed bond-lengths. I believe you used LAMMPS for most of your MD simulations. Have you started using GROMACS? I know you could learn it quickly, but my concern is the short deadline. Also, GROMACS requires tabulated values for the Mie potential, which are a lot slower than the native LJ function.

6d- We might combine options a) and d) to at least include an example for when cut-offs are problematic. Your ITIC results for Option d) might not be compatible since I do not believe we have the trajectories/coordinates that we need for MBAR. But they could be used to check the GROMACS energies and pressures.

I could be persuaded in favor of either Option c) or d). Option c) is a bit more "plug and chug" but would certainly meet the FPE expectations and would help for IFPSC. Option d) is more interesting to me but is also less straight-forward, i.e. we could run into some unforeseen obstacles. Option d) could certainly be a publication for a different journal though. Let me know which direction you are leaning.

Rich

[ramess101]

Rich,

6c- Regarding your first concern, I think we can count on our machines in Akron. we have a powerful Tesla K80 gpu machine that we want to put in good use. We also have 3 GTX Nvidia gpu machines as well as two 24-core and 16-core cpu machines. What kind of computational resources do you have at NIST? Regarding your second concern, I agree that there will be a learning curve, since I never used GROMACS to get useful data, however I have tested GROMACS and I have a few scripts that facilitates running GROMACS simulations. I wouldn't worry much about this issue. I remember somebody in our video conference mentioned that GROMACS is hard to modify. How hard is that? It would be very beneficial to add Mie and other potential models to GROMACS. Tabulated potentials cannot be fast and accurate at the same time. We might need to look into modifying GROMACS, at least for validating our results.

6d- It's a good idea to combine a and d. Yes, we can only use my ITIC results to check GROMACS P and U since they do not have large trajectory files, only a few snapshots.

In general, Dr. Elliott and I are leaning towards option d. It has more impact on our long-term goals and it's an opportunity for me to get my hands dirty with MBAR.

S. Mostafa Razavi

[ramess101]

Mostafa,

6c- I have access to 28 cores on a newer cpu machine and access to 16 cores on an older cpu machine. The speed is not the main issue it is the number of cores. For EMD viscosity runs we typically run 30-50 replicates at each state point. This is why I think it might be difficult to simulate a large number of compounds with fewer than a couple hundred cores. 6c- I agree that adding other potentials could be very beneficial. I believe someone added a patch that allows for Mie potentials in GROMACS, but I think it essentially converts your Mie input into a tabulated potential, since it did not appear to be any faster.

6d- I am also leaning towards this option, mainly because it is more interesting, feeds our long-term project, and could have a high impact. My main concern is whether we will have something ready by August. Then again, I think that this approach is novel enough that it does not necessarily need to go in the FPE special issue.

Let me think about this some more over the weekend. We should also talk about these options with Jeffrey Potoff in our next video conference.

Rich