Thermodynamic Integration with MACE in LAMMPS

LAMMPS
,
1

Hello,

I was reading a paper (doi.org/10.1021/acs.jpclett.5c02980) that lays out a method for computing solvation free energies using MACE potentials and thermodynamic integration. The paper used OPENMM, but I wanted to see if it is something that is achievable in LAMMPS. The novelty of the paper was to use the thermodynamic integration variable lambda to affect the neighbor list by artificially scaling the distances between the solvent and solute. In the simplest form this looks like:

dij’ = dij + lambda * rcut, with i being an atom of the solvent and j being an atom of the solute, other pair interactions are unaffected.

So at lambda = 0, dij’=dij, and at lambda = 1 dij’ = dij+rcut > rcut, therefore the solute and solvent don’t see each other.

I suspect this isn’t possible in LAMMPS, an alternative approach I thought of was to use TI to move the solute from solution (position A) to the vacuum (position B). The energies and forces would just be an interpolation between the 2 states i.e. U = (1-lambda)UA + lambdaUB ect. with the solvent maintaining the same micro state for the solute in both cases. I don’t know if I explained that clearly enough, but in any case it involves tying lambda to the energy or 2 separate energy/force evaluation.

I’d welcome any feedback on if these 2 methods are possible in LAMMPS or if there is an alternative I have not considered. Thank you.

2

@mlacount please note that MACE is not part of LAMMPS and developed and maintained separately, any issues related to MACE itself you have to discuss with them.

LAMMPS supports a whole lot of free energy methods, specifically through the COLVARS and PLUMED packages, e.g steered MD, ABF, Metadynamics and more), many of which should be applicable to compute solvation free energies rather independent from the model used for interactions. I don’t know anything about specific requirements or restrictions of MACE, so that is also a question for the MACE developers.

3

Thank you for your reply. If I understood correctly this implies that at least method 1 I described would need a modification of MACE. I thought this was probably the case, but wanted to check my assumption.

My second method may still be viable, however the way MACE works provides a complication. MACE only gives you the total energy and the net forces on each atom, there is no way to decompose the forces into solvent-solute and solvent-solvent for example. Except if you run the calculations twice with 2 different geometries both with and without the solute.

This is the feature I don’t know if LAMMPS (or COLVARS/PLUMMED) supports. Basically to compute the energies and forces twice with the same pair style but different set of atoms active (say 1 to 99 then again with 2 to 100), and have the result of mixing them, U = (1-lambda)*UA + (lambda)*UB (similar with the forces on each atom), used in the molecular dynamics/thermodynamic integration. My first thought on how to do this would be to use either pair_style hybrid/overlay or pair_style hybrid/scaled, both with the same pair_style but different pair_coeff to select out the solute in one case. However, I couldn’t quite get it to work.

The manual on pair_style mliap says “A single pair_coeff command is also required. The first 2 arguments must be * * so as to span all LAMMPS atom types.” and that seems to a hardcoded requirement so 2* 2* is not valid.

4

What you are proposing makes physically no sense, sorry. But it also looks to me like you need to make a more thorough survey of what methods are available to compute (solvation) free energies.

Here are your options in LAMMPS to look into and that I can remember (there may be more):

  • there is thermodynamic integration using compute ti. This would require that you are using a pair style that is compatible with with fix adapt and specifically you would have to have support for a “scale” parameter. The fx adapt documentation has a list of supported pair styles. changing pair style mliap to apply a scale parameter and adding to an extract() function is not very difficult. there are many precedents.
  • you can model alchemical transformations with fix alchemy (using two partitions!). this is completely independent from any modifications or special features of the pair style, it requires a careful setup of the two partitions, though.
  • you cannot partition the pair style computation (that is unphysical with a manybody potential and would result in garbage data), but you can scale resulting forces and energies for a single, full calculation using pair style hybrid scaled in combination with an atom style variable as scaling factor, i.e. you can have a different scaling factor for each atom, e.g. for different atom types.
  • you can use a completely different free energy method like steered MD or ABF or metadynamics where you instead of “growing” or “vanishing” a molecules you pull it across a liquid/vaccum interface. these are supported by both COLVARS and PLUMED.
5

Thank you for your reply. To clarify I am trying to do a technique I’ve seen in literature (https://doi.org/10.1063/1.4983754) to compute the solvation free energy via a thermodynamic cycle: create a cavity, insert the solute, collapse the cavity. The first and third step are fairly straightforward, its the middle step I am trying to find a solution for.

  • As you said for the first option it would require a pair style that supports a scale parameter, which could be done but would require a change of the code. I’ll keep this in mind as an option.

  • This seems to be the closest to what I am looking for, except that according to the documentation: “simulation starts with the same atom positions across both partitions”. So trying to use it to move the solute spatially is not an option. It had occurred to me to use dummy atoms that alchemically transform into the solute, but my understanding is that dummy atoms can’t be used with the pair_style mliap.

  • I think I can see how this could work, with a slight modification to the code. I would need to modify lammps to accept arbitrary ‘pair_coeff’ rather than just ‘pair_coeff * *’.

  • This was my original approach, I was just running into difficulties trying to decouple the perturbation in the interface with the free energy. In short, it caused a jump in free energy in the direction of the pull, so depending on whether I was pulling the solute into or out of solution I’d see a strong hysteresis in my free energy curves. This led me to explore alternatives like alchemy/thermodynamic integration, to avoid having to deal with the interface, but I could return to trying to solve the hysteresis problem instead.

6

I already explained that this is wrong and will not work. This condition is enforced for a good reason.

It looks to me that you need some in-person assistance from somebody that has more patience than I have and and more experience than you have.