Disagreement between Cassandra and Towhee

[ramess101]

Here are the results from simulations performed on both Cassandra and Towhee. The LJ results assume the same sigma and epsilon as the ab initio potential (this is mainly just to provide some level of confidence):

In Cassandra we use the exact Fortran code for the two-body potential, while in Towhee we use a tabulated potential.

These results are obtained using a 1.4 nm cut-off WITHOUT tail corrections. The choice to not include tail-corrections is because Towhee does not allow for tail-corrections with a tabulated potential.

All simulations were performed using 2800 molecules, where the initial configuration has 2400 in the liquid phase and 400 in the vapor phase.

Towhee simulations use an equilibration of 20,000 MC cycles and production of 100,000 MC cycles. Cassandra simulations use an equilibration of 2,000,000 MC moves and production of 2,000,000 MC moves (2,000,000/2800 is approximately 714 MC cycles, so the Cassandra simulations are considerably shorter. This might be the cause for the deviation, but the plots appear to have converged despite the relatively short equilibration period.

Towhee results correspond to 8 replicates whereas Cassandra has only 3 replicates.

[ramess101]

We have already validated that Cassandra and Towhee compute the same energies for the same initial configuration. So the energies should not be the issue.

[ramess101]

To give an idea as to how close the Towhee and Cassandra non-bonded energies are:

[ramess101]

I am now running the Lennard-Jones potential in Cassandra to see if I can capture the correct VLCC curve. If not, I will know that there is something wrong with my version of Cassandra or how I am implementing GEMC. The move probabilities might be worth looking at more closely.

[ramess101]

The Cassandra and Towhee results for the LJ potential (using sigma and epsilon that correspond to the Helium ab initio potential) provide indistinguishable results (especially considering these LJ results are from a single, relatively short, replicate simulation). Note that I have Towhee LJ results using both the native LJ and the tabulated (truncated with 1.8 nm cut-off). The results are the same, suggesting that Towhee's interpolation scheme is not to blame.

So this does not provide any clear evidence as to why our Towhee and Cassandra results are so different for the ab initio potential.

[ramess101]

I ran some additional simulations at higher temperatures in Cassandra, and the critical point appears to be much higher than in Towhee:

[ramess101]

The ratio between the Cassandra and Towhee liquid densities is not constant with respect to temperature:

7 K : 1.04333032031 8 K :1.06769895416 9 K :1.09348579397 10 K :1.12277809576 11 K :1.16582397014

Suggesting that this is not some molecular weight or units issue in the final output.

[ramess101]

@Kaiveria

I found a small discrepancy in the Towhee and Cassandra energies. My python script that is used to generate the tabulated potential in Towhee had a less precise conversion factor for bohr to Angstrom. Specifically, the value I used was 0.529177 while in potentials.f90 for Cassandra we used 0.529177249. I noticed this a long time ago, but I never really considered it a potential source for error. However, the deviation between the energies is slightly larger than I would have expected. It is still small enough that I would not anticipate this to impact the GEMC results, but the systematic bias (albeit small) could accumulate.

Here are the percent deviations, notice that the "Different Conversion Factor" curve has a -0.0003% bias over the entire attractive region. This means that the Towhee simulations are slightly less attractive than they should be (a negative percent deviation in the attractive region, where U is negative, means that Towhee minus Cassandra is positive). A less attractive potential would, in theory, lead to a narrower coexistence curve (the hard sphere, zero attractive potential, has no coexistence). However, I still strongly doubt it would cause this large of a deviation.

Here the absolute percent deviations are plotted since they do not have a singularity as r = sigma (U=0) and they show that we have reached machine single precision:

So I am going to rerun the Cassandra simulations using the wrong conversion factor (I know this seems backwards, but it should be faster than running Towhee with the correct conversion factor).

[ramess101]

This small deviation in energy does not seem to be the cause of the VLE deviation. My Cassandra simulation results with the wrong conversion factor are indistinguishable from those with the correct conversion factor.

[ramess101]

It should be possible to determine whether the Towhee or Cassandra results are correct based on the LJ results. Notice that Towhee has a narrower vapor-liquid coexistence curve (VLCC) while Cassandra has a wider VLCC than that of the LJ fluid.

Here is what the ab initio and LJ potentials look like:

So the LJ potential has a more attractive shoulder than the ab initio potential and is much harder at close distances. Other than that they are quite similar.

Since the ab initio potential is slightly less attractive at moderate distances and more repulsive at close distances, I would expect that the VLCC curve would be narrower than that of the LJ. This would suggest that Towhee is more correct. However, this type of analysis is not infallible as there are other subtle differences between the potentials.

[ramess101]

David Kofke was able to verify that the Towhee results are more accurate than those from Cassandra:

Andrew ran some simulations last night, and his findings follow below. He did two runs using a fast but approximate potential, and two using the full pair potential.

In short, it looks like Towhee is giving correct results.

Dave

run 1: T=10K, N=200+200 atoms, approximate potential, rho_liquid=rho_vapor=25 mol/L, cutoff 15A liquid density: 62.25(11) mol/L = 249.00(44) kg/m^3 vapor density: 2.75(5) mol/L = 11.0(2) kg/m^3 liquid pressure: -0.106(8) MPa vapor pressure: 0.192(3) MPa

run 2: T=10K, N=1000+1000 atoms, approximate potential, rho_liquid=rho_vapor=25 mol/L liquid density: 63.20(10) mol/L = 252.8(4) kg/m^3 vapor density: 2.90(4) mol/L = 11.60(16) kg/m^3 liquid pressure: 0.21(3) MPa vapor pressure: 0.1992(17) MPa (negative liquid pressure is a finite size effect, other properties are not super sensitive)

run 3: T=10K, N=300(liquid)+100(vapor) atoms, full (2012) potential, rho_liquid=63, rho_vapor=2.75, cutoff 10A liquid density: 62.28(12) mol/L = 249.12(48) kg/m^3 vapor density: 2.74(4) mol/L = 10.96(16) kg/m^3 liquid pressure: -0.07(3) MPa vapor pressure: 0.194(2) MPa (approximate potential is indistinguishable from full potential)

For good measure, one more run completed this morning:

run 4: The full potential with N=1000+1000, liquid density: 63.04(9) mol/L (252.3(4) kg/m^3) vapor density: 2.69(6) mol/L (10.8(2) kg/m^3) liquid pressure: 0.19(4) MPa vapor pressure: 0.188(3) MPa

The liquid density appears to match Towhee (254), certainly nowhere near Cassandra (287). Vapor density is harder to discern, but Towhee looks like 11 kg/m^3 and Cassandra looks like 5 kg/m^3, so Towhee appears to win for vapor as well.

[ramess101]

Things I have tested recently that did not fix the problem:

1) I verified that the linear interpolation works well in Towhee. In fact, it appears to be much more reliable than the cubicspline, i.e., the energy for an initial configuration was in closer agreement with the true potential.

2) I used an external LJ function in the same way that we use the external ab initio potential. But the results were the same as using the internal LJ.

3) I tweaked the frag1.dat file values.

4) I changed the potential from Mie to LJ.

5) I moved where V and V_diff are called to inside some if statements

6) I compared the NVT energies for Towhee and Cassandra. They are not markedly different.

7) I played around with the precision in the energy computation.

[ramess101]

We found the bug! Again, it was a units problem. The energy should be output in atomic units, not Kelvin. The LJ potential converts epsilon to atomic units (amu A^2/ps^2) prior to running the computation. This is a factor of 0.831 or 1.2 different. Essentially, we were running our simulations at epsilon*1.2 because we were not multiplying our energies by 0.831. This explains why our Tc was too high.

With this epsilon we would get this LJ curve:

So the Cassandra results are to the left of the LJ curve (for a 1.2 * 11 epsilon).

[ramess101]

The Cassandra jobs have completed their equilibration stage and have begun the production stage. The uncertainties are still relatively large, but this figure shows the good agreement with Towhee achieved:

[ramess101]

This plot is useful for visualizing the vapor phase:

[ramess101]

Here are the final plots after allowing the Cassandra runs to finish (although they are still considerably shorter than Towhee they are long enough for these preliminary purposes, i.e. 1400 without tail corrections):