ramess101 / helium_ab_initio Goto Github PK

This study by Siepmann's group demonstrates a robust approach for estimating the critical point:

The question is, do we need this type of precision? For example:

1. Do we need to simulate that close to the critical point
2. If so, we need to use density distribution functions. If not, we probably don't need to use density distributions.
3. How important are the MC moves?
4. How important is the liquid to vapor ratio?
5. Do we need to use very large systems?
6. Do we need to vary the cut-off?
7. Do we need to use their approach for estimating uncertainty?

I have written a python script that determines the required LJ epsilon and sigma to achieve the actual tail corrections for the full potential:

`

rcut = 14 #[A]

U_tail_rcut = U_tail(rcut)
P_tail_rcut = P_tail(rcut)

print('Cutoff distance: '+str(rcut)+' [A]')

print('Actual tail corrections, U and P:')
print(U_tail_rcut,P_tail_rcut)

eps_LJ = 11. #[K]
sig_LJ = 2.64 #[A]

ULJ_tail = lambda eps,sig,rcut: 4.*eps*(1./9. * sig**12. / rcut**9. - 1./3. * sig**6. / rcut**3.)
PLJ_tail = lambda eps,sig,rcut: -8.*eps*(2./3. * sig**12. / rcut**9. - sig**6. / rcut**3.)

print('Initial guess LJ tail corrections, U and P:')
print(ULJ_tail(eps_LJ,sig_LJ,rcut),PLJ_tail(eps_LJ,sig_LJ,rcut))

Utol = 1e-10
Ptol = 1e-10

for i in range(1000):

    sigma6 = -rcut**3./(8.*eps_LJ)*(P_tail_rcut + 12.*U_tail_rcut)
    sigma12 = sigma6**2 
    sig_LJ = sigma6**(1./6.)
    eps_LJ = U_tail_rcut/(4./9.*(sigma12/rcut**9.) - (4./3.)*(sigma6/rcut**3))
    if i%10 == 0:
        devU = ( ULJ_tail(eps_LJ,sig_LJ,rcut) - U_tail_rcut ) / U_tail_rcut *100.
        devP =  ( PLJ_tail(eps_LJ,sig_LJ,rcut) - P_tail_rcut ) / P_tail_rcut *100.
                
        if devU < Utol and devP < Ptol: 
            print('Number of iterations required and converged epsilon [K] and sigma [A]')
            print(i,eps_LJ,sig_LJ)
            break

print('Final LJ tail corrections, U and P')
print(ULJ_tail(eps_LJ,sig_LJ,rcut),PLJ_tail(eps_LJ,sig_LJ,rcut))

assert devU < Utol, 'Not precise enough energy tail correction'
assert devP < Ptol, 'Not precise enough pressure tail correction'`

This is the output for rcut=14 Angstrom:

Cutoff distance: 14 [A]
Actual tail corrections, U and P:
(-1.3594884274740304, 7.7702233663389269)
Initial guess LJ tail corrections, U and P:
(-1.8095240801733516, 10.856981755313758)
Number of iterations required and converged epsilon [K] and sigma [A]
(690, 0.0028662102721911643, 10.037020621578602)
Final LJ tail corrections, U and P
(-1.3594884274740289, 7.7702233663453999)

Notice that epsilon is 0.0029 and sigma is 10.

Here it is for rcut = 18:

Cutoff distance: 18 [A]
Actual tail corrections, U and P:
(-0.7054775255750203, 3.8236498317381247)
Initial guess LJ tail corrections, U and P:
(-0.8514046675265038, 5.108411055694046)
Number of iterations required and converged epsilon [K] and sigma [A]
(370, 0.00037622243257559964, 14.421128645581341)
Final LJ tail corrections, U and P
(-0.70547752557501986, 3.8236498317398473)

Epsilon is now 0.00038 and sigma is 14.4. This is interesting but also as expected since the ab initio potential converges to zero more quickly than the LJ.

The three-body module (he3body.f90) is now working properly. Here is a validation run:

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.