Free energy differences due to potential scaling

[g-bauer]

There are various free energy methods that utilize scaling of intermolecular interactions. For example one can compute the chemical potential by inserting an additional particle/molecule into the system. This insertion can be difficult if the molecule is big or the system is dense. A way to circumvent these difficulties is to slowly grow the molecule by coupling the intermolecular parameters such that the molecule is an ideal gas molecule at the starting state and a fully interacting molecule in the end state.

The free energy can then be computed e.g. using thermodynamic integration. It would be nice if Lumol would support scaling of potential functions as well as computation of properties that are usually of interest when computing free energies, such as the derivative of the energy w.r.t. the coupling parameter.

The easiest way to couple the potential is a linear coupling (simply multiplying the energy with a coupling parameter) but that might not always be ideal. It would be nice to allow for implementation of coupling methods as well as some support to conduct free energy computations in the input file.

[Luthaf]

I agree that this functionality is something we want to have, because if is widely useful: one can use it to implement free energy computations, grand canonical MD, fractional component MC, ...

How do you think we can implement this in the code ? We need a way to store the scaling parameter that allow us to use it when computing energy and to update the scaling parameter during the simulation.

One possibility would be to store the parameter inside a shared, reference counted pointer, that would be shared between a potential and the part of the simulation responsible for the update (MC move, MD propagator/control algorithm, ...).

[g-bauer]

Could another trait on top of Potential work? Given a scaling parameter, one could define the energy function (or fall back to linear scaling by calling to Potential energy), force function and the derivative w.r.t the scaling parameter (since that's the property usually used for FEP computations). Not sure where to store this maybe in a struct similar to Interactions? And also not sure where coupling parameter(s) is/are stored. Often, one needs more than a single coupling parameter, i.e. if coulomb and vdw contribution are separately coupled.

[Luthaf]

Yep, another trait could work too, storing the scaling parameters in Interactions. This would be the best way to do it if you need the code quickly.

However, I would like as much as possible to reuse existing code, instead of adding specific code path for every new features. If we can not implement this feature within the current framework, maybe a change in the framework will allow us to implement both the current version and this new feature? If not, then we can add new code path :smile:.

The other solutions that I can think of using the current framework are:

(1) Sharing the parameter by using a Rc pointer. This way, some code can update it (the MC Move) and some other code can use it (the energy computation code):

struct ScaledPotential<T> {
    scaling: Rc<f64>,
    potential: T
}

// Implement energy/force/...

// Add function for the derivative w.r.t. the scaling parameter.

// also add the parameter in a MC move
struct ScaleMCMove {
    scaling: Rc<f64>,
}

impl ScaleMCMove {
    fn new(potential: ScaledPotential) -> Self {
        // The same parameter is shared between the MC move and the potential
        // When the move updates the paramater, the potential changes
        ScaleMCMove {
            scaling: potential.scaling.clone()
        }
    }
}

(2) Using the restriction code for this, by adding a Scaled case:

enum PairRestriction {
    // Currently
    None,
    IntraMolecular,
    InterMolecular,
    Exclude12,
    Exclude13,
    Exclude14,
    Scale14(f64),
    // Add this
    Scaled {
        scaling: f64,
        particle: ParticleKind,
    }
}

// Use this to implement scaling whatever interaction
// the particle have with any other particle in
// PairRestriction::information

// The the code can dynamically change the scaling
// paramater value by accessing the PairInteractions

(1) has the advantage of inserting cleanly in the current code, and the inconvenient of needing more work from the user: the user needs to add both a potential to the system and a Move to the simulation sharing the same Rc.

It is also harder to make it work with coulomb interactions, but I think the code could just dynamically adjust the charge of the particle(s) to scale the Coulombic energy?

(2) Might need some changes from the current state of the code, and the way we compute the excluded part of RestrictionInfo (we'll need more that the bond distance to compute it), but feels a bit cleaner to me.

We also get Coulombic interactions for free here, as coulomb already work with PairRestriction.

The main downside I see is that we can not scale intra-molecular interactions, and that we are effectivley limited to pair and coulomb potential. Would this be a downside for you?

---

Concerning the derivative w.r.t. the saling parameter, it should be possible to use forces here, right? Like for linear scaling, the derivative is - potential.force().

[g-bauer]

Concerning the derivative w.r.t. the saling parameter, it should be possible to use forces here, right? Like for linear scaling, the derivative is - potential.force().

For linear scaling, the derivative w.r.t. the scaling parameter (for the energy) would simply reuse potential.energy() (since the scaling function is simply "scaling_factor times energy", but that only works for linear scaling. For LJ potentials e.g. soft-core scaling is very common and in such cases we need a special routine for energy, force, derivative, ...

I think it will be difficult to build the utilities into the current framework if we want it to be as general as possible (which is to be discussed). I have no experience about which utilities are offered by other codes such as LAMMPS and DLPOLY but in Gromacs (at least from the API point of view) there are utilities for van der Waals, coulomb, restraint and bond (also mass and temperature) perturbations which can be done independently, i.e. there is a "scaling parameter" for each of them.

I am not sure how to implement this in the current framework. My first guess would be to create scaled versions inside Interactions as you said in your post. Not sure about restrictions since cases may overlap.

[Luthaf]

OK, if we want/need to also change the bonds parameters and the mass and so on, the solution with restriction is not very good.

I think the proposition (1) could still be used: simply sharing the scaling parameter between all the possible users of it (interactions, MC moves, ...) This would also work with mass/charge, as the MC move could update them on the fly. Temperature seems harder.

[g-bauer]

Sure. Personally I only use(d) vdw and coulomb perturbations and have no experience with the other stuff. vdw + coulomb (intermolecular) would be a good starting point, I think.