So, there are a few things to address here.

(1) When you say 'Gibbs free energy change', this implies a starting state and an ending state for the system. What is the end state that you are interested in? Often people are interested in taking a molecule out of its surroundings and putting it into gas phase, where it has negligible interaction energy with its surroundings.

In the case of removing a molecule to gas phase, the LJ-potential is an approximation of dH, though not necessarily a good one if the pressure-volume term is significant (i.e. LJ potential does not tell you how much work was done to expand the molecule into gas phase). Even if you assume the PV term is negligible, you still need to estimate -TdS to get dG.

Unfortunately, entropy change usually cannot be measured with the LJ calculation of a single snapshot because a single snapshot does not tell you what other conformations are accessible to the molecule. It could be the case that there are dozens of other conformations with similar or better energies (in which case the entropy of the current state is favorable). Or it could be the case that all other conformations have very high energies (in which case the entropy of the current state is unfavorable). Entropy of a molecule can be thought of as measuring how many 'choices' the molecule has.

(2) Are you running a molecular dynamics or Monte Carlo simulation? What sort of molecules are you simulating? Entropy can be estimated from simulations because they give you an idea of how many other conformations are accessible.

In some cases you can use experimental data to estimate the entropic cost of your molecule's state change. If you tell us more about your setup, perhaps we can advise you further.

Best,

Hassan

According to what I have learn from my life long experience is that  the people in reality are able to calculate some how energy plus the entropy, and put them together  to read the Helmholtz characteristic function,  F= U-TS,  then   they argue that the stored strain energy density  in solids is  negligible small compared to all other energy sources especially electronic in nature.  That means  G  is about equal to F.  Actually, In the statistical Thermodynamics books by Fowler and Guggenheim openly and sincerely They admit that the calculations of Enthalpy from the first principle is almost a forbidden task for solid.  One has the same problem in the variational  formulation of surfaces and interfaces in deformed solids using irreversible thermodynamics of isothermal changes as advocated by solid Ogurtani unless one assumes that  the subject  matter solid is exposed to small deformations acting as an linear elastic media obeys the Hook's law or it is hyperelastic for the  finite deformations as introduced  by Late Professor Truesdell.   Actually the main reason is our ignorance to find a Lagrangian function to  represent  the deformation, which  involves  in addition to space variables but also the stress tensor and derivatives of stress tensor, similar to the Lagrangian function Representing Helmholtz characteristic function (See: Eshelby in Solid State Physics Series V.5, 1957).

Thanks for your answers.

I will try the suggestions of Hassan and get back to you!

Check the Widom particle insertion algorithm. or, use the grand canonical ensemble to "estimate" it! A good starting point may be the book by Frenkel and Smit.

All statistical mechanics and thermodynamics formulations rely on the Energy  U  and Helmholtz free energy F characteristic functions  , respectively, if they  deal  with the infinitesimal adiabatic and isothermal processes taking place in close systems. That creates very serious problem  if one deals with the deformed solids exposed to external surface tractions and body forces. As far as the entropic part  is concerned  according to Callen ; for the most solids the effect of stress on the entropy is almost negligible .

However, in one of my published work I have shown how this stress  effect on the Gibbs and Helmholtz free energy densities   can be formulated using the Latent Heat Tensor referred to the unstrained  state. 

dear Fernando Vallejos-Burgos,

A first suggestion is to make variable temperature calculations of the potential, U. Plot U(T) vs the inverse temperature (1/T). Presumably this will be shaped somewhat linear. The slope is related to dH and the intercept to dS. I assume this is quite a routine but utilized frequently in experiments.

The second suggestion is to use thermodynamic integration. This one is approached strictly with molecular simulations. The method is well described in Book Understanding molecular simulation : from algorithms to appl...

In brief, use a decoupling parameter, l, for the potential expression.

That is U(l)=(1-l)U, where l in [0..1]. For every l, compute U in NPT at a few nanoseconds time scale (depends on the system). Make a histogram of the average, U, vs l. Integrate this to get the free energy change.

Unfortunately there is no direct simple thermodynamic connection between the slope of U(T) versus 1/T plot and the enthalpy of the system as suggested by Dr. Gotzias, which can be easily realized if one looks at the fundamental differential equations of thermodynamics as called by Gibbs (SEE: Guggenheim, p.28: Eqs: 1.33.1- 1.33.4). Actually the partial derivative of U(T;V) with respect to T under the constant pressure results the specific heat capacity under constant pressure Cv. Namely: d (U(T,V) /d T)V = (TdS /dT)V = Cv

On the other hand dU(T,V)/d(1/T)V = -T2 d[U(T,V)/dT]V =- T2 CV

As I have mentioned before one should define the boundary conditions or constraints of the process in terms of [P or V] and [ T or S]. For a isothermal isobaric process we choose {T, P} where the characteristic function to be selected should be Gibbs free energy G(T, P, composition). In the case of isovolumic process (isochoric) the chose would be F(T, V) for the isothermal case, which is mostly employed in statistical thermodynamics because F= U -TS where one has to compute U (T:V) and S(T:V) for a given volume, Temperature and configuration of the particles involved. Here mostly one mostly skips temperature (T=0 ) rather tries to fine the most probable configuration of the particles which optimize the Internal energy of the system rather then the Helmholtz Free energy of the system.

For the adiabatic system one should distribute the give energy over the particles in terms of transitional kinetics energies randomly in order to thermalized the system.