Reviewer 2 Comment 1

[ramess101]

@mrshirts @jpotoff @msoroush

section 3.4: this case study ignores a paper where a Mie force field (also lamda=16) was earlier proposed for cyclohexane, with (I guess) essentially the same results. Differences should appear only because the objective function is slightly different here. The conclusion of this work is that the here proposed force field is „the most accurate“. Omitting the previous study is (at least) awkward, because the authors seem to be aware of the line of developments from these authors, as the reference list shows. Is there a reason for omitting the study mentioned above from the list of „most reliable force fields from literature. Ref 22,24,27,55,56,57“ ?

The reviewer never clarifies what paper has a Mie force field for cyclohexane. @jpotoff @msoroush Any ideas which paper the reviewer is referring to? Who are "these authors" that are included in our reference list? I can't find one for TAMie (the only other Mie family we cite). I know there are parameters obtained with SAFT-gamma that use fractional lambda values and some single-cite cyclohexane Mie models. But I don't think we have any data for those models.

[ramess101]

@jpotoff @msoroush

OK, I am pretty confident that the reviewer is referring to TAMie. I didn't notice that there were some TAMie parameters reported for cyclohexane in the aldehydes and ketones paper. They did not reparameterize cyclohexane in the "individualized" study so I assumed they didn't have any parameters.

Their parameters are quite similar to ours. I wish I would have known this beforehand, it sort of steals our thunder.

TAMie: epsilon = 69.568 sigma = 3.8967 lambda = 16

MiPPE: epsilon = 69.7 sigma = 3.902 lambda = 16

Of course, it is still impressive that our first iteration (starting with TraPPE) yields similar parameters to the TAMie optimal:

First iteration: epsilon = 70 sigma = 3.89 lambda = 16

I could include that TAMie point on our heat map to emphasize this point.

TAMie reported tabulated property values too, so I can certainly include those in our deviation plots.

[jpotoff]

@ramess101 That's unfortunate that we missed their cyclohexane model. I should have known better as I've seen that paper before, but my brain was a bit foggy last semester. Still, the paper is about a method and an optimization process, so this isn't a deal breaker. We include the TAMie data, and use it as validation of the method. We start from TraPPE, optimize parameters, and get to TAMie (about). I'm wondering what our data say about the deviation between the predictions of the TAMie parameters and experiment. Is it the same is what is published, or do we see slightly different results? That might be a good comparison to include in supporting information if they match. If they don't...Maybe still throw it in supporting information. If we don't have it in the heat map, it would only take a day or two to run it (and I think we should).

BTW, I do have a concern about this being an apples to apples comparison, though. If we subjected the exp-6 and/or the TAMie potential to the same thoroughness of calculation (20 independent runs), I would expect to see some shifts in the data.

[ramess101]

@jpotoff

Still, the paper is about a method and an optimization process, so this isn't a deal breaker. We start from TraPPE, optimize parameters, and get to TAMie (about).

I agree. It certainly doesn't minimize the influence of our methods, but the culmination with MiPPE cyclohexane parameters losses a bit of its punch. I need to modify the story we are telling a bit in the Case Study section too.

I'm wondering what our data say about the deviation between the predictions of the TAMie parameters and experiment. Is it the same is what is published, or do we see slightly different results? That might be a good comparison to include in supporting information if they match. If they don't...Maybe still throw it in supporting information. If we don't have it in the heat map, it would only take a day or two to run it (and I think we should).

I will let you know once I have these results. Even if our results differ slightly from the literature TAMie values, this could be attributed to finite-size effects (considering their alternative DeltaHv approach) and partially the difference in torsion potential. But I am pretty confident they will agree within the combined uncertainties.

BTW, I do have a concern about this being an apples to apples comparison, though. If we subjected the exp-6 and/or the TAMie potential to the same thoroughness of calculation (20 independent runs), I would expect to see some shifts in the data.

The main reason of including Exp-6 and TAMie results was to meet the JCED requirement that we compare our results with preexisting literature values.

[mrshirts]

The main reason of including Exp-6 and TAMie results was to meet the JCED requirement that we compare our results with preexisting literature values.

Right, so hopefully we can do this, and show that the method works . . .

BTW, I do have a concern about this being an apples to apples comparison, though. If we subjected the exp-6 and/or the TAMie potential to the same thoroughness of calculation (20 independent runs), I would expect to see some shifts in the data.

And that could be pointed out.

[ramess101]

@mrshirts @jpotoff @msoroush

I have not finished modifying the manuscript, but this is what I am planning:

[ramess101]

@mrshirts @jpotoff @msoroush

This is what the deviation plots look like with TAMie included (grey circles):

[ramess101]

@mrshirts @jpotoff @msoroush

I am going to include these figures in the SI. Clearly the MiPPE (and to a lesser extent the first iteration of MiPPE) results agree very nicely with the TAMie reported values. Interesting that rhol is the largest deviation. The purported finite-size effects for DeltaHv at low T are nonexistent.

[ramess101]

@jpotoff @msoroush

I'm wondering what our data say about the deviation between the predictions of the TAMie parameters and experiment. Is it the same is what is published, or do we see slightly different results? That might be a good comparison to include in supporting information if they match. If they don't...Maybe still throw it in supporting information. If we don't have it in the heat map, it would only take a day or two to run it (and I think we should).

OK, after doing some digging we did compute phase properties for the TAMie parameters in our refined second iteration heat map scan, at least if we round the parameters to within 3 digits of precision (which is not uncommon). Specifically, TAMie reports epsilon of 69.568 K and sigma of 3.8967 Angstrom while we predicted values for 69.5 or 69.6 K and 3.895 or 3.900 Angstrom. The values shown here as "TAMie, This work" are for 69.6 K and 3.9 Angstrom.

I think this is good enough to show that our results are consistent. I could run for the EXACT same epsilon and sigma, but I think we are well within the noise at this point. Also, I don't have access to my system at NIST where I ran these calculations originally.

We could always run a full set of GCMC simulations for TAMie (20 replicates, etc.). But I think that is overkill. Any objections?

[ramess101]

@jpotoff @msoroush

But if you and Mohammad want to run 20 replicates with (exact) TAMie, I will certainly include those results instead of mine. My GCMC-MBAR results are admittedly not the best values to report because they are estimated from a single set of simulations performed with the first iteration theta parameter set. Mohammad would probably need to play around with the best mu values for TAMie. And he might want to run with 3.5 nm (?)

[msoroush]

@ramess101 I can run the TAMie by tomorrow. What is the exact parameter you want me to run? And what box size?

[ramess101]

@msoroush

OK, perfect! Let's use the exact values that TAMie reported:

Weidler et al. reported using box sizes between 40,000 A^3 (for propanal) and 81,000 A^3 (for octanal), so I would say that we should use the 35 Angstrom box length instead of 30 Angstrom box length (since this is only 27000 A^3).

[msoroush]

@ramess101 I also used their reported torsion

[ramess101]

@msoroush

Good point. So you have their modified torsion potential coded in GOMC?

[msoroush]

@ramess101 I used the cos series but I converted it using C0 = C0 - (C1 + C2 + C3)

[msoroush]

@ramess101 using this cosine series Utorsion = C0 + C1 (1+cos(phi + delta1)) + C2 (1+cos(2*phi + delta2)) + C3 (1+cos(3*phi + delta3)) I used this value:

C_n          n          delta
-5339        0           ---
6840         1           0.00
3509         2           0.00
63           3           0.00

[ramess101]

@msoroush

Perfect. Interesting that they (both TAMie and TraPPE) didn't just correct the C0 term in their table instead of making it look like it is a completely different equation.

[msoroush]

@ramess101 have you tried to plot the torsion energy for the equation and parameter used in TAMie? Is the minimum correct?

[msoroush]

@ramess101 @jpotoff C1 term must be negative. Otherwise we face the same problem. I am sampling boat configuration and no crown! You might put back your old statement about wrong torsion value.

[ramess101]

@msoroush

I have not, I just assumed that because they changed the equation that it fixed the problem. I will go ahead and put back in the old statement unless you get this to work.

So should we simulate TAMie with the same torsion that we used for MiPPE?

[msoroush]

This is TAMie torsion with original parameter:

This is TAMie torsion with modified parameter (-C1 term):

This is alkane torsion parameter:

[msoroush]

@ramess101 @jpotoff

Is that okay if I use modified torsion of TAMie?

[ramess101]

@msoroush @jpotoff

I think whatever you used for MiPPE should be ok for TAMie. The TAMie equation and parameters are supposed to be the same as TraPPE. And so if changing the equation still requires making c1 negative, I think it would be simpler to just use what we did for MiPPE previously (use alkane).

[msoroush]

@ramess101 @jpotoff Okay. I will use alkane torsion parameters.

[ramess101]

@mrshirts @jpotoff @msoroush

Included this statement that MiPPE and TAMie are very similar. Also wanted to make it clear that we are not purporting that MiPPE is better than TAMie since TAMie used a slightly different objective function (relative squared error for only rhol and Psat) and experimental data set (DIPPR):

[ramess101]

@msoroush

Do you expect the TAMie results by tomorrow? I think that is the last thing that I need before I can resubmit. Thanks

[msoroush]

@ramess101 I am still working on fixing the chemical potential. I will have results for one replica. I will set it to run for 20 replica, by tonight but not sure when they will be finished. When is the deadline to resubmit it?

[ramess101]

@msoroush

That's fine. The deadline is next Monday but I just didn't want to cut it too close, to give us time in case the results don't come out as we expect.

[msoroush]

@ramess101 Do you want the single replica results? Wrong torsion parameters messed up my schedule :D

[ramess101]

@msoroush

Haha, those torsions are going to be the death of us someday...

Sure, pass the single replica results along when you have a chance. That way I can make sure they look reasonable compared to Weidler et al.

[msoroush]

@ramess101 here is the data for single replica. Simulations are not finished, I had to kill them to start other simulations. So, the quality might not be as good as we want, specially at lower temperature. vle.zip

[ramess101]

@msoroush @jpotoff

I am a little bit concerned about our TAMie replicates. I know this was just a single (short) replicate, but the deviation in Psat and rhov is way greater than what Weidler reports.

[ramess101]

@jpotoff @msoroush

OK, I looked into this a little more. And I think your trends for TAMie are actually correct. The figures that I plotted previously with my own GCMC-MBAR predicted curves for TAMie were not for the correct sigma. I plotted results for eps = 69.6 K (which is the closest I have to TAMie of 69.568), but I was using MiPPE sigma (0.3902) instead of TAMie sigma (0.3895 or 0.3900 are the two closest I have to 0.38967).

Here are my results for sigma of 0.3895:

Here are my results for sigma of 0.3900:

If you mentally interpolate the Psat and rhov errors for the two sigma values, you will see that my TAMie results actually do agree with your TAMie results. In other words, the Weidler et al. Psat and rhovap values are quite different. Could this be related to the torsions? I think we will need to anticipate our TAMie results not matching Weidler et al. and provide some explanation for this in either the review or the manuscript.

[ramess101]

@jpotoff @msoroush

The TraPPE results provide additional evidence that the torsions might be to blame for this discrepancy. Our Psat results agree with those of Yiannourokou (who also used the alkane torsion potential) but disagree with Keasler (who purportedly used the TraPPE cyclic potential). We don't really agree with either for rhovap, but I think this might explain the few percent shift in Psat for TAMie.

[ramess101]

@msoroush @jpotoff

So after these 20 replicates finish, I think we might want to try simulating TAMie using the TraPPE cyclic potential (but with a negative c1 value).

[msoroush]

@ramess101 for MIPPE data, are we using 30 or 35 A data?

[ramess101]

@msoroush

Good point. These MiPPE values were for 30 just because that is what we used in the optimization. But I could use 35 instead to be consistent with our TAMie values. The results do not differ by several percent though (except at high T) so that doesn't really explain the discrepancy.

[msoroush]

@ramess101 I am looking at my TraPPE results and I see the same discrepancy in vapor density. I think, the best way to sort it out is to run NPT simulation in gas phase using MD code, for each torsion model (alkane torsion and TraPPE torsion) and compare the predicted densities.

I am not sure how with wrong torsion parameter, they can get these results. In GOMC if I used the wrong torsion (TraPPE torsion), my liquid density will be much higher than what they reported.

[ramess101]

@msoroush

Yeah, I don't know how they could get their results with the reported torsions. But maybe if it really is a sign error for C1, the results would shift slightly?

[msoroush]

@ramess101 Here are the results for alkane torsion and modified torsion (-C1 term). I dont see any changes in the data.

cyclohexane.tar.gz

[ramess101]

@msoroush

Strange. So if it is not the torsions, then maybe there is some subtle difference in the methods, codes, etc.

What about the cutoff length? Are we using 10 A or 14 A? TAMie and TraPPE both use 14 A cutoff.

[msoroush]

@ramess101 @jpotoff

The problem is that vapor density results for Keasler does not match with Yiannourakou, either.

Using rcut 10 vs 14 would not affect the results as long as we use Long Range Correction. We can look at the branched alkane results and compare it with TraPPE results.

[ramess101]

@msoroush @jpotoff

True, but Keasler used GEMC and didn't report uncertainties, so I am more interested about the TAMie comparison. I expect the TAMie Weidler results are quite reliable and were computed with a similar GCMC approach (although they did use a different method for computing pressure as mentioned in Comment #29).

I agree that with long-range corrections you typically get indistinguishable values for rcut of 10 and 14, but in some cases the assumption of g=1 at rcut=10 may not be valid (rings?).

[ramess101]

@msoroush @jpotoff

My vote would be that we repeat the single replicate (at least) with rcut of 14 for TAMie, just to be safe. The box length is 35 A, so rcut = 14 should not be a problem there.

[msoroush]

@msoroush @jpotoff

My vote would be that we repeat the single replicate (at least) with rcut of 14 for TAMie, just to be safe. The box length is 35 A, so rcut = 14 should not be a problem there.

Okay, I will set it up.

[ramess101]

@msoroush

Thanks! Sorry for being so demanding... I just want to make sure we leave no stone unturned

[msoroush]

@ramess101 Just to be clear, you want me to repeat the cyclohexane simulation for TAMie model using rcut 14 A and 35 A cubic box?

[ramess101]

@msoroush

Yeah. And using the alkane torsion is fine too.

[msoroush]

@ramess101 Do you have time to setup GROMACS to run NPT simulation for liquid and vapor phase?

[ramess101]

@msoroush

Sorry, I actually do not. I have zero computational resources (not set-up yet at NREL and no access at NIST). So I am kind of useless and I definitely can't get anything done by Monday. 👎

[msoroush]

@ramess101 You ran the initial simulation for TraPPE in GOMC? Would you please plot the densities predicted by GOMC, Keasler, and Yiannourakou?