Dear lammps users,
I’m willing to use LAMMPS to compute the Henry constant of heptane in kerogen (a polymer like material).
As a test case, I first want to validate my input script on a simpler case, by considering pure heptane at 700 K. The heptane molecules are represented with the OPLS-AA force fields. For the time being, I do not consider Coulombic interactions (all charges are set to 0).
Thus, my system is the following : I have a pure heptane fluid and I’m willing to compute the excess chemical potential (also called solvation free energy) of a tagged heptane molecule with the other (untagged) heptane molecules. Then, using the Henry law, I could determine the Henry constant from the excess chemical potential.
To determine the excess chemical potential of the tagged molecule, I choose a Finite Difference Thermodynamic Integration (FDTI) scheme where d_U / D_lambda is computed with the compute/fep module.
In a prior step, using lj/cut pair style for unbounded interactions, I got an excess chemical potential of ~1.30 kCal/mol, which is close to the excess chemical potential obtained from a Peng-Robinson equation of state (1.40 kCal/mol).
Due to the singularity when lambda = 0, I’d rather use of softcore potential. However, when switching to a lj/cut/soft potential, within the same conditions, I now get an excess chemical potential of 0 kCal/mol. These results have been obtained with both the stable version of LAMMPS from the 15 may 2015 and the last version of lammps (8 sept 2015).
Here is the part of my input script related to the definition of the pair coefficients and the compute/fep routine :
variable TK equal 700
# Scaling variables
variable lambda equal 0.1
variable dlambda equal 0.002
variable s_cc equal 3.5 # sigma C-C for CH3
variable s_hh equal 2.5 # sigma C-C for CH2
variable s_ch equal sqrt({s_cc}*{s_hh}) # sigma H-H for Halk
variable e_cc equal 0.066 # epsilon C-C for CH3
variable e_hh equal 0.03 # epsilon C-C for CH2
variable e_ch equal sqrt({e_cc}*{e_hh}) # epsilon H-H for Halk