Reviewer 3, Comment 6: uDep from single molecule

[mostafa-razavi]

Eq. 7 and 8 are written as indefinite integrals, but they not useful in this form. Eq. 8 is used to construct Eq. 10, where it becomes a definite integral. It would be simpler to write them as definite integrals from the start.

I would suggest that the derivation would be simpler and less confusing if the derivation targets full free energy

A(rho,T^IT) = A_id(rho,T^IT) + integral[rho'=0,rho (d(A-A_ig)/drho) drho']
beta A(rho,T) = beta A(rho,T^IT) + integral[beta'=beta^IT,beta (d betaA/dbeta) dbeta']

Splitting the free energy into ideal and non-ideal parts is necessary for density integration (because the full integrand would diverge at rho=0), but it's unclear what it accomplishes for the temperature integration, except perhaps to make the results incorrect: the authors later (in section 4.1) describe subtracting out the intramolecular energy (taking U_ig = U_bond + U_intra).

This seems wrong. The U_bond and U_intra contributions in a simulation at some finite density are not the energy an ideal gas would experience. The ideal gas contribution might instead be determined from a simulation of a single molecule where U_ig = U_tot = U_bond + U_intra.

Stated another way, the path from coexisting vapor to coexisting liquid might be represented as

(ρV,Tsat) ==1==> (0,Tsat) ==2==> (0,TIT) ==3==> (ρL,TIT) ==4==> (ρL,Tsat)

Rather than explicitly accounting for the change in ideal gas free energy in stage 2, the authors instead attempt to avoid it entirely by subtracting it out in stage 4. But the energy used for this purpose seems not to be the ideal gas energy (as would have been seen in stage 2), but the intramolecular (including bonding) energy seen in stage 4. I don't see how that can be right. If the authors are computing ideal gas energy for stage 4 (using a single molecule in a box), the they should state this explicitly. And, in practice, the whole thing would be simpler if the ideal gas contribution to the full free energy were included along (in stage 2) instead of trying to subtract the same from stage 4.

[mostafa-razavi]

Our response:
The reviewer raises a point that has worried us for several years. For starters, it is not entirely correct to say that the zero density limit corresponds to the ideal gas energy when considering a chain molecule like nC12. An ideal gas is composed of point masses, without any non-bonded intramolecular energy, whereas the nC12 that we simulate does include non-bonded intramolecular energy at zero density. That low density non-bonded intramolecular energy is implicitly accounted for through the virial coefficients, although we never actually compute the non-bonded intramolecular energy directly. We simply compute B2(T) and solve for the Helmholtz energy that is implicit in the temperature dependence of B2. We can surmise, of course, that the non-bonded intramolecular energy is non-zero at zero density because at least a few conformations must lead to wrap-around contacts that are non-zero. (Incidentally, we have performed a fairly detailed analysis of this effect in the context of discontinuous molecular dynamics and thermodynamic perturbation theory for step potentials. That work was part of an AIChE presentation several years ago. The effect seems interesting to us, but frankly the audience did not seem interested at the time, so we set it aside.) The question of whether to subtract the intramolecular non-bonded energy is most easily addressed by comparing results from ITIC to GEMC. The figure below shows a comparison of vapor pressures and liquid densities computed with and without the non-bonded intramolecular energy. Once again, we compare deviations from REFPROP values as a baseline to make the magnitudes of the discrepancies more clear. Clearly, failure to subtract the non-bonded intramolecular energy leads to results that are inconsistent with GEMC. Briefly, we conclude that GEMC leads to equilibration of the intermolecular interactions between molecules, without interference from intramolecular interactions, neither bonded nor non-bonded. It is an esoteric point but quite apparent from the ITIC analysis, which should serve to make the ITIC analysis that much more interesting as a complement to the other alternatives. We have incorporated these figures into the supporting information and refer to them in the context of subtracting the non-bonded intramolecular energy as a necessary step in the ITIC process.

[ramess101]

@mostafa-razavi

Make sure that the SI has a clear caption for this figure to clarify the symbols. You might even want to include it in the review.

[ramess101]

@mostafa-razavi

Should we just remove the Ungere rhov data? I suspect there is actually a units issue as they are off by 1000%.

[ramess101]

@mostafa-razavi @jrelliottoh

Is this correct?

An ideal gas is composed of point masses, without any non-bonded intramolecular energy, whereas the nC12 that we simulate does include non-bonded intramolecular energy at zero density.

I thought an ideal gas also includes non-bonded intramolecular interactions. I don't see why the bond strecthing, angle bending, and torsional interactions would be include without the non-bonded 1-4, 1-5, etc. interactions. Just because we refer to these interactions in force field lingo as "non-bonded" doesn't really make them any different than the standard bonding terms. For example, ab initio calculations of a single molecule yield ideal gas properties while implicitly including non-bonded interactions. Right?

[mostafa-razavi]

Should we just remove the Ungere rhov data? I suspect there is actually a units issue as they are off by 1000%.

I already did.

TraPPE-dodecane

Mie-dodecane

TraPPE-isohexane

[ramess101]

I already did.

OK, but did you also check to see if you just needed to divide their rhov by 1000?

[mostafa-razavi]

@ramess101 @jrelliottoh

Some observations and conclusion from the above figures:

1. The difference between uDep calculated using the two different approaches is on average around 1.7 %.
2. uDep from single molecule simulations is more negative, therefore the Psat is more positive.
3. As seen in the above figures, this difference causes a significant deviation in Psat and rhoV
4. Psat deviation increases by decreasing saturation temperature. Psat difference between the two approach is 10-15 % for Tr=0.45 and 1-4 % for Tr=0.85.
5. I think, the deviation between the two approaches is worse for C12, because C12 is bigger than isohexane.
6. Interestingly, using single molecule approach solves the Psat inconsistency between GCMC and ITIC for Mie-C12 completely, while rhoV has similar deviation in the opposite side.
7. If we trust the single molecule approach I think we should recommend in the text that for large molecules, the single molecule approach is preferred.
8. We should also replot Figure 10, 11 (only C12 and iC6 results), and Figure 12 based on the single molecule approach.

What do you think?

[ramess101]

Just a quick note about GCMC. I have been a little concerned and confused why Potoff uses the nonbonded energy rather than total energy when reweighting histograms. Is there a chance that the GCMC results actually have a small bias because of this?

…

On Friday, December 21, 2018, mostafa-razavi ***@***.***> wrote: @ramess101 <https://github.com/ramess101> @jrelliottoh <https://github.com/jrelliottoh> Some observations and conclusion from the above figures: 1. The difference between uDep calculated using the two different approaches is on average around 1.7 %. 2. uDep from single molecule simulations is more negative, therefore the Psat is more positive. 3. As seen in the above figures, this difference causes a significant deviation in Psat and rhoV 4. Psat deviation increases by decreasing saturation temperature. Psat difference between the two approach is 10-15 % for Tr=0.45 and 1-4 % for Tr=0.85. 5. I think, the deviation between the two approaches is worse for C12, because C12 is bigger than isohexane. 6. Interestingly, using single molecule approach solves the Psat inconsistency between GCMC and ITIC for Mie-C12 completely, while rhoV has similar deviation in the opposite side. 7. If we trust the single molecule approach I think we should recommend in the text that for large molecules, the single molecule approach is preferred. 8. We should also replot Figure 10, 11 (only C12 and iC6 results), and Figure 12 based on the single molecule approach. What do you think? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#26 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AWUvhIVEVFS7dby9KLDmlpG75kkWo3y_ks5u7KUbgaJpZM4Zc2YE> .

[ramess101]

@mostafa-razavi

Can you send me the reference for Ungerer TraPPE dodecane? I can't find it. I just want to verify the units for rhov.

[mostafa-razavi]

Optimization of the anisotropic united atoms intermolecular potential for n-alkanes

[mostafa-razavi]

112-40-3.txt

[ramess101]

@mostafa-razavi @jrelliottoh

I am working on revising this response so that it is easier for the editor and reviewer to see that we addressed their concerns appropriately.

[ramess101]

@mostafa-razavi @jrelliottoh

Mostafa and I think that we can include the intra method for now for a few reasons. 1) It is simpler 2) Single molecule simulations are problematic in MD 3) It works for small systems 4) Including this discussion provides some insight so that the user will know why they need to use single molecule for large molecules.

[ramess101]

@mostafa-razavi @jrelliottoh

We are planning on updating the C12 results in the figures such that they use single molecule and mentioning in results that this was necessary for a larger flexible molecule. Same with the Figure 12 system.

[mostafa-razavi]

Here is the new Figure 10. C2 results are updated

[mostafa-razavi]

I think we should remove TraPPE-C12-Gromacs because it's still using Uintra method.

[ramess101]

@mostafa-razavi

Looks pretty awesome to me. Just make sure to mention how we use single molecule for Cassandra. I think Gromacs C12 is still fine to include without correction. Note that a single molecule MD simulation typically requires some stochastic sampler. So I don't even think we would necessarily recommend MD for the single molecule. We could instead just use the Cassandra single molecule to correct Gromacs as well, if desired.

[ramess101]

@mostafa-razavi

How hard would it be to correct Gromacs with Cassandra single molecule? I think we are justified here because we could just say that single molecules were always run with MC instead of MD.

[mostafa-razavi]

How hard would it be to correct Gromacs with Cassandra single molecule? I think we are justified here because we could just say that single molecules were always run with MC instead of MD.

It needs a little digging. I'll try that if I find some time.

[mostafa-razavi]

I think we should use single molecule approach for isohexane as well. There is ~10 % difference in Psat fro Tr=0.45. here is the plot using single molecule isohexane

[mostafa-razavi]

@ramess101 Any comment on the above figure?

[ramess101]

@mostafa-razavi

Yeah that is fine. We don't really know which is better at Tr=0.45, but single molecule should be better.

[ramess101]

@mostafa-razavi

Wordsmithing:

Udep is computed with Eq (30) for n-dodecane and isohexane by performing single molecule simulations with Cassandra, while Udep is computed with Eq (29) for other compounds.

[mostafa-razavi]

Udep is computed with Eq (30) for n-dodecane and isohexane by performing single molecule simulations with Cassandra, while Udep is computed with Eq (29) for other compounds.

It's not entirely true, because we use GOMC for Mie-C12. Maybe I forgot to mention that.

How about
Udep is computed with Eq (30) for n-dodecane and isohexane by performing single molecule simulations with MC, while Udep is computed with Eq (29) for other compounds.

[ramess101]

@mostafa-razavi

If it was different for each system, I wouldn't even say "with MC."

But so are some of these ITIC results (not just single molecule simulations) with Cassandra and some GOMC? Or are they all GOMC?

[mostafa-razavi]

They are all Cassandra except for Mie-C12 which is GOMC and filled black symbols which are GMX.

[ramess101]

@mostafa-razavi

OK, just don't say anything in caption about how you computed single molecule simulations. Is it clear in the manuscript or at least the SI which simulation packages were used for the higher density simulations? That is important information to include.

[ramess101]

@mostafa-razavi

Do you have the isobutane single molecule plots so that I can include them in my response?

[mostafa-razavi]

Do you have the isobutane single molecule plots so that I can include them in my response?

[mostafa-razavi]

Is it clear in the manuscript or at least the SI which simulation packages were used for the higher density simulations? That is important information to include.

It's mentioned in SI.

[ramess101]

@mostafa-razavi

Small typo in Udep discussion (molecules misspelled as "moleucles"):

[ramess101]

@mostafa-razavi

Can you upload the most recent .tex file? I would like to work on section 4.1 some and see if you like how it comes out.

[ramess101]

@mostafa-razavi

Also, the current version should say, "An underlying assumption of Eq. (29)" not Eq 10. Eq 10 just says U-Uig, but how we compute Uig is the question.

[mostafa-razavi]

Can you upload the most recent .tex file? I would like to work on section 4.1 some and see if you like how it comes out.

Updated in repository.

[ramess101]

@mostafa-razavi

OK, I ended up just using an old version since I was rewriting most of it anyways. Here is what it looks like right now. Any comments? If you like it I will send you the .tex file.

[mostafa-razavi]

@ramess101 It's perfect. Thanks. Please send me this tex file so I include it

[mostafa-razavi]

@ramess101 I'll fix the broken links

[mostafa-razavi]

@ramess101 Are you also gonna modify the response accordingly?

[ramess101]

@mostafa-razavi

OK, the links might not be broken when you copy it to your file since I used the old names.

[ramess101]

@mostafa-razavi

Are you also gonna modify the response accordingly?

Yes.

[ramess101]

@mostafa-razavi

We need to reference this paper for the discussion of MC or SD instead of MD:

Merz PT, Shirts MR (2018) Testing for physical validity in molecular simulations. PLoS ONE 13(9): e0202764. https://doi.org/10.1371/journal.pone.0202764

[ramess101]

@mostafa-razavi

I made a few more modifications to help it read better. I also provided a rough estimate for the percent error in Psat if intra is used with large molecules. Let me know if this statement is not accurate.

Here is the tex file (remember to add the reference above):

intramolecular.zip

[ramess101]

@mostafa-razavi @jrelliottoh

Here is my revised response to the reviewer. Let me know if something looks wrong.

Revised_Udep.docx

[ramess101]

@mostafa-razavi

Here is the bibtex for that additional reference of stochastic dynamics instead of MD:

@Article{10.1371/journal.pone.0202764,
author = {Merz, Pascal T. AND Shirts, Michael R.},
journal = {PLOS ONE},
publisher = {Public Library of Science},
title = {Testing for physical validity in molecular simulations},
year = {2018},
month = {09},
volume = {13},
url = {https://doi.org/10.1371/journal.pone.0202764},
pages = {1-22},
number = {9},
doi = {10.1371/journal.pone.0202764}
}

[mostafa-razavi]

@ramess101 Got it. thanks

[ramess101]

@mostafa-razavi

Please post a picture of the revised section once you have recompiled the PDF so I can verify all the links are correct.

[mostafa-razavi]

[ramess101]

@mostafa-razavi

Looks good