Excess chemical potential of heptane with compute FEP

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

Normally the lj/cut and lj/cut/soft should give the same answers when lambda = 1. I will test this in detail ASAP.
Meanwhile, maybe you could try to use lambda as the parameter being adapted between your initial and final states (from 0 to 1 I suppose) and keep the epsilon and sigma fixed (that’s the idea behind the soft pair styles).

Agilio

2

Dear Agilio,

First of all, thank you for your help.

Normally the lj/cut and lj/cut/soft should give the same answers when lambda = 1. I will test this in detail ASAP.

As you mentioned the two potentials lead to the same value of U_1 - U_0 when lambda = 1.

Meanwhile, maybe you could try to use lambda as the parameter being adapted between your initial and final states (from 0 to 1 I suppose) and keep the epsilon and sigma fixed (that’s the idea behind the soft pair styles).

I also did some tests using the fix adapt/fep, but it gaves me the same results. I have reworked a little bit my data to compare the results obtained from FDTI with both Lennard-Jones and soft-core potential. I provide you the spreadsheet in attached.

It is organized as follow :

  • In the left hand side of the spreadsheet are reported different thermodynamic outputs (kinetic, potential and total energy, fep1 fep2 and fep3 as defined in the lammps manual), for both soft-core and Lennard-Jones potential. The next columns corresponds to the weight coefficient for the integration (wi) and du/dlambda obtained from fep1 and fep2 ;
  • Then, on to right hand side of the sheet are reported the results : excess chemical potential, and for the problem I’m interested in, the fugacity in both lammps and « user friendly » units ;

revue-fep-fdti-ti.ods (32.7 KB)

3

Hi Julien,

I’ve looked at the spreadsheet you sent with the results and have two questions:

  • 700 K may be a stretch for OPLS-AA, which is parameterised at milder conditions. Could it be that the divergence when using lj/cut when lambda -> 0 leads to a value closer to experiment just by accident? The curve using the soft-core potential looks much better behaved.
  • how did you change your parameters when using lj/cut? Linearly epsilon leaving sigma constant?

Best,

Agilio

4

Hi Agilio,

  • 700 K may be a stretch for OPLS-AA, which is parameterised at milder conditions. Could it be that the divergence when using lj/cut when lambda -> 0 leads to a value closer to experiment just by accident? The curve using the soft-core potential looks much better behaved.

What was struggling me with the results for the softcore potential was that the factor exp(mu_res / kt) was close to unity, which (if the results are correct) would indicate that the fluid is close to an ideal gas. But, I will relaunch new simulations to check that point. Similarly to the examples you provided in the lammps distribution, I will try to determine the solvation free energy of heptane in liquid water at 300K.

  • how did you change your parameters when using lj/cut? Linearly epsilon leaving sigma constant?

Yes, this is exactly what I did : varying epsilon at constant sigma. In that case, lambda corresponds to the prefactor used to scale the epsilon, lambda varying from 0 to 1. As it can be seen from the spreadsheet, I made a denser sampling between lambda = 0 and lambda = 0.1 due to the strong variation of the slope of the curve in that range.

Thanks for your help,
Best regards.