Argon

[ramess101]

Because I am having problems with the absolute pressures for hexane, I am now going to try with argon since this compound is used in the workshop: https://github.com/GOMC-WSU/Workshop/blob/master/GCMC/argon. This should facilitate debugging as I can compare my intercept value with the value they reported.

Here are the histogram plots from the seven state points that I simulated with GOMC:

[ramess101]

First, I verify that I get the same rhov and rhol points as reported by Mick et al.:

[ramess101]

Here are the results from my analysis compared with the reported intercept value from the workshop file

https://github.com/GOMC-WSU/Workshop/blob/master/GCMC/argon/patch/phinput.idat

Slope for ideal gas is 1, actual slope is: 0.999996090033 Intercept for absolute pressure is:-2.43235462657

@msoroush Can you verify for me that the value of -2.4061 in the workshop argon phinput.idat is the correct intercept? If it is, I need to figure out why my intercept is so different.

[msoroush]

@ramess101, Here is my results. You can try the hist file that I attached and use the chemical potential listed in pvt.idat to calculate the molecule number and ln of partition function.

argon.tar.gz

[ramess101]

@msoroush

Thank you for getting back to me. My intercept actually only differs from yours by about 0.02. This would not lead the large deviations in absolute pressure that I am observing. Furthermore, if I use the histogram that you provided, I get the exact same intercept:

Slope for ideal gas is 1, actual slope is: 0.999997682922 Intercept for absolute pressure is:-2.45489349921

So this demonstrates that the lnZ0 term is not the reason my pressures are so off. This was very helpful, as I can now look somewhere else for the bug in my code.

[ramess101]

@msoroush

I continue to have an error in my pressures. To figure out what the problem is, it would be very helpful if I had your Zliq values that are used to compute the pressures.

I tried to compute Zliq myself using your code. Specifically, I followed the README instructions in the workshop (https://github.com/GOMC-WSU/Workshop/tree/master/GCMC/argon) to run the patch/program files on my own histograms, but I kept getting this error when I tried to run ./pvt.out:

PHASE VERSION 8.0 PVT MODE At line 178 of file pvt.f90 (unit = 2, file = 'input_fsp2.dat') Fortran runtime error: Bad real number in item 5 of list input

We could try to figure out why I am getting this error, but it would probably be easier and faster if you just ran your code on my histograms (see attachement). Could you analyze these histograms and provide me with the Zliq and pressure values?

Thanks

ram9his.zip

[msoroush]

@ramess101 try to delete weights.dat and input_fsp2.dat files. Then try again following README file. I had the same issue, this would solve your problem.

[msoroush]

@ramess101 here are the results for pressure, liquid and vapor densities and Z. pcov generate the pressure in tp.dat and trho.dat for densities, but it does not calculate the Z. You can double check my calculation for Z. argon.dat.zip

[ramess101]

@msoroush

Thanks. Your code seems to be working for me now. I verified that our weights are the same for each of the 7 simulation points, i.e. my code gives the same values as the weights.dat file. But for some reason there is a nearly constant deviation of -0.69 in our lnZliq values.

[ramess101]

@msoroush

I was able to find the lnZliq values in phase.dat

[ramess101]

@msoroush

Sorry, I think I confused you when I kept saying Z. I was referring to the partition function (specifically, lnZliq), not the compressibility factor.

[msoroush]

:D

[mrshirts]

-0.69 = -ln 1.9937 . . . ? Loose factor of 2 anywhere?

[ramess101]

Hmmm... so maybe somewhere I am missing a -ln(2)

[ramess101]

@mrshirts

OK, all I can think of is that I need to divide by 2 somewhere to account for the liquid and vapor phases. Currently this is how I get the lnZliq:

Deltaf_ij = mbar.getFreeEnergyDifferences[0] lnZliq = Deltaf_ij[nTsim:,0] #nTsim is the number of simulated temperatures. nTsim: yields the Deltaf_ij values for only the VLE (non-simulated) mu,T points

press = kT/V*(lnZliq - lnZ0)

But we get the correct pressure if we use:

press = kT/V*(lnZliq - ln(2) - lnZ0)

Since @msoroush and I got the same lnZ0 values for the ideal gas intercept, this suggests that:

lnZliq_correct = ln(Zliq/2)

Since Zliq = exp(Deltaf_ij), this would mean that the exp(Deltaf_ij[nTsim:]) needs to be divided by 2. @mrshirts can you think of why this would be true? Again, my only guess is that somehow this needs to account for the two phases.

These equations from Issue #1 might help

Note that Ci is Deltaf_ij, so exp(Ci) is exp(Deltaf_ij), which is equal to the sum over all E and N of p.

[ramess101]

@msoroush @mrshirts

Also, note that the -ln(2) term not only fixed my argon pressures but also my hexane pressures. So I think this really was the problem.

Thanks again for helping me resolve this issue! Now we just need to be able to explain it.

[mrshirts]

It's not entirely clear to me where this difference is occuring from. If @msoroush ISN'T seeing this, and is using the same histograms, can you compare their code? Is it possible the error is in that code?

[ramess101]

This snippet from https://github.com/GOMC-WSU/Workshop/blob/master/GCMC/argon/patch/program/phase.f90

    `if(n(1,ifile)%flength(i).lt.nmid(tnum)+split)` then
        Zgas = Zgas + exp(-ylog)

                   nbelow = nbelow + exp(-ylog)
                    do icomp = 1,ncomp
                       ngas(icomp) = ngas(icomp) + &
           n(icomp,ifile)%flength(i)*exp(-ylog)
                    enddo
                    energy_gas = energy_gas + &
                         e(ifile)%flength2(i)*exp(-ylog)
    else 
        Zliq = Zliq + exp(-ylog)

                    do icomp = 1,ncomp
                       nliq(icomp) = nliq(icomp) +&
           n(icomp,ifile)%flength(i)*exp(-ylog)
                    enddo
                    energy_liq = energy_liq + &
                         e(ifile)%flength2(i)*exp(-ylog)
            end if`

clearly shows that Zliq and Zgas are only added to when N is greater and less than Nmid, respectively.

[ramess101]

If @msoroush ISN'T seeing this, and is using the same histograms, can you compare their code? Is it possible the error is in that code?

@mrshirts

I am very confident in their code. My original pressures had some serious issues that were clear red flags that my code had the problem. For example, my logP vs 1/T trend was very strange. Also, at low Tsat where the vapor is nearly ideal (because vapor pressure is very low) my pressures were not converging to the ideal gas value.

Note that these results were using my MBAR-GCMC code. I have not tried to fix my histogram reweighting code because it had some other issues that were not worth debugging.

[mrshirts]

Are there the same number of molecules in the liquid and vapor phase? Do you get a different correction factor if the number is uneven?

[ramess101]

@mrshirts

The number of molecules in each phase is quite different because the volumes are equal.

[mrshirts]

I can't think of a physical justification for "there are two phases, so divide by 2". What is the ratio of the volumes that are used?

[ramess101]

@mrshirts

The volumes are the same.

My explanation is the following:

For traditional histogram reweighting, in Equation 15 we compute exp(Ci), which is the same as exp(Deltaf_ij) in MBAR lingo, by summing the probabilities of all the bins.

As shown in @msoroush code above, this summation is broken into two terms, Zliq and Zgas, based on the Nmid (i.e., low N adds to Zgas and high N adds to Zliq). So I think essentially what I have been using is Zliq+Zgas, which at VLE equilibrium these are the same, and thus we have exp(Deltaf_ij) = 2*Zliq.

[mrshirts]

Ah, OK. If the simulation uses the total energy and total N over all boxes (i.e. semigrand), then yes, it's the total partition function. The partition functions of the phases aren't equal at equilibrium, but pressures should be, so for NmuT, then P V = /propto Z, so for equal volumes, then you would have the total Z = 2*liquid Z. You could check this by running a simulation with a different (fixed) volume for the gas, and reconstructing what the correction term is in that case.

[ramess101]

@mrshirts

If the simulation uses the total energy and total N over all boxes (i.e. semigrand), then yes, it's the total partition function. The partition functions of the phases aren't equal at equilibrium, but pressures should be, so for NmuT, then P V = /propto Z, so for equal volumes, then you would have the total Z = 2*liquid Z.

OK, that makes more sense.

You could check this by running a simulation with a different (fixed) volume for the gas, and reconstructing what the correction term is in that case.

That would be interesting. I think we can revisit this idea later, once we have everything else working well.