## Feature #1652

Feature #1332: Supporting multiple end states instead of just A and B

### Decide how to represent multiple lambda states internally

**Description**

One of the obstacles to solve is what the internal representation of the different ends states should be, and how one described the intermediate states that interpolate between these states. This thread breaks out that subdiscussion.

### History

#### #1 Updated by Michael Shirts over 5 years ago

Some questions that will need to be answered. It will require an analysis of what the goals of the multiple state code are (see the parent discussion)

- Are there multiple lambda vectors, controlling where one is with respect to each lambda state?
- If so, how do these lambdas get defined in the .mdp?
- Are these from the same root state (A->B, A->C, A->D), consecutive (A->B, B->C, C->D), or pairwise (all N choose 2 pathways between N endstates?)
- Is it possible to be halfway from A to B and halfway from B to C? i.e. do N intermediates define an N-dimensional space, and one can be anywhere in that space, or can the system only be along edges?

#### #2 Updated by Thomas Ullmann over 5 years ago

We have two aims which motivate the design of the of the implementation that I'm working at. Each of these aims includes a number of features that we would like to have. I tried to also incorporate the features that Michael has identified as important here and in a previous Email.

### A) We want to consider systems with multiple (potentially many) variable sites that can interconvert between different forms during a simulation, e.g., titratable sites that can uptake or release protons.¶

- The sites should be allowed to have any number of forms to model protonation forms, tautomers, redox forms, binding of other (small) ligands, slowly interconverting conformations ... and any combination of them at the same time.
- Closely related to the previous point, we want to have the possibility to consider not only protons but any number of other small ligand species, e.g., electrons, metal ions, entire small-molecule ligands. Each site form has than an entry for each ligand type that specifies how many ligands of this type are bound.
- Again closely related, membrane proteins are often situated between two compartments with different electrochemical potentials for a given ligand type, e.g., the proton. In presence of an electrostatic transmembrane potential, there is also an energy contribution to a binding reaction that arises from the removal of a charged ligand from the respective compartment to which the site has access and its placement in the receptor environment. Therefore each site form needs to have a record that specifies to which compartment the site has access.
- It should be possible to use different variants of lambda dynamics / interconverting between the site forms (further discussed below).
- It should be possible to let different types of force field terms interconvert/switch separately, e.g., Coulomb and Lennard-Jones interactions.

The theory to treat the complex thermodynamics described by A.1.-3. is fully developed and described partly in a previous paper on the simulation code GMCT [[http://dx.doi.org/10.1002/jcc.22919]]. If someone is interested in the details there is also a detailed, yet rather simple derivation in my PhD thesis.

The current state of the implementation of my "lambda site" module already has a mechanism to realize A.4. Allowing for A.5. also is, I think, straightforward and noted on my ToDo list.

An brief explanation of the terminology used in my implementation maybe necessary here to make sure that the rest of my post is understandable. Each variable site/"lambda site" comprises a site-specific number of site forms. Each site form has an associated lambda variable that measures its contribution to the mixture of forms. A form is fully present if its lambda value is equal to 1, and fully absent if its lambda value is equal to 0. The sum over the lambda values of all site forms

should ideally be close to one. A physical state is signified by lambda = 1 for a single form.

### B) On top of A) we want to perform free energy calculations on these systems where the system is forced along a transformation coordinate between two (or more?) end states.¶

- states (end states and intermediates between them/alchemical states) of the transformations/reactions considered in these calculations can be quite complicated.

a) part of the variable sites may be restricted to subsets of their forms

b) the system may be restricted to a range of (collective) coordinate values

c) a) and b) are defined separately for all end states that contribute to

an alchemical intermediate - it should be possible to consider whole sets of different transformations/reactions

with a user defined free energy calculation method and target statistical error, that is

ideally without the need for postprocessing the simulation results

I think that B should codewise be treated separately from A.

Consequently, there would be a lambda variable on top of the lambda variables of the sites that controls the transformation. For clarity, I will distinguish this lambda from the lambda values of the site forms by naming it "w" for weight.

In the simplest case there your free energy calculation considers the free energy difference between two forms, 0 and 1, of a single site, w would be given by

w = lambda_0 = 1 -lambda_1. In the general case, there is no such one-to-one correspondence.

More than two states, as defined above, could also easily be considered at the same time (each simulated separately with each global lambda/w having a fixed value for a given state).

I'm not sure whether one would want to evolve also the global lambda values that control multiple lambda sites

as continuous variables. Currently, It seems to me that this would not be necessary. Otherwise, formalism and implementation may become quite complicated.

A working proof of concept for these ideas is implemented in another simulation code for models with discrete states (GMCT).

Are there multiple lambda vectors, controlling where one is with respect to each lambda state?

Yes, provided I understood the question correctly. Each site will have such a vector, possibly with separate lambda vectors for different force field terms to address A.5.

The global lambda/w value(s) would of course also be tracked.

Are these from the same root state (A->B, A->C, A->D), consecutive (A->B, B->C, C->D), or pairwise (all N choose 2 pathways between N endstates?)

In my formalism (and similarly in previous, more specialized versions), the energies of the different site forms are defined relative to one freely chosen reference form. This does, however, not hinder free energy calculation schemes under B) that consider multiple end states at once, sequences of transformations between pairs of states or anything else that one could want.

Is it possible to be halfway from A to B and halfway from B to C? i.e. do N intermediates define an N-dimensional space, and one can be anywhere in that space, or can the system only be along edges?

This leads to a discussion on different lambda dynamics variants. Several of them have been published, two more are under investigation by Helmut Grubmueller, Gerrit Groenhof, Serena Donnini, Plamen Dobrev and myself.

A bit on my variant:In my view, this variant has several advantages, where the first two directly address the above question:

- the system can directly access each form from each other form. This does however not mean that one will not be able to do free energy calculations that interconvert between sequences of pairs of the forms.
- The system will be allowed to have contributions from more than two forms at the same time to profit as much as possible from the sampling enhancement through lambda dynamics. An adjustable bias potential will give the user control over the number of forms that can contribute to the mixture at the same time.
- any number of forms can be considered
- the sum of lambda values of the forms of a site is guaranteed to be a constant, typically equal to 1.

Helmut, Gerrit and Serena and Plamen explore a variant that evolves the system along linear edges of grids spanned between different forms. As far as I can see it, this variant will become very complicated if one has more than four site forms, with considerable complexity already for more than two forms. In addition, it requires additional constraints that ensure that the weigths of the forms sum to 1. It may however have advantages in certain situations, and my variant may have weaknesses that I did not yet realize.

A modular design of our lambda site implementation will allow us to explore the different variants (and possibly others), without much coding effort. The invariables of lambda dynamics are covered by a common lambda site object (calculation of the energy terms of each site form required to compute the forces along the lambda coordinates, setting up the lambda site from the input topology ...). Aspects that may need to be treated by different variants, like how the site evolves along the lambda coordinates, or which MD integrator to use are exchangeable and referenced by the lambda site object via pointers. The modularity requires the variable parts to have as few dependencies as possible. For example, the MD integrator should be a function object without any internal state variables which just accepts input coordinates and forces and returns the updated coordinates without knowing about the detailed nature of the coordinates (I think this has been proposed also by others on the developers list for MD integrators in general, If I remember correctly by David van der Spoel).

The architecture of the lambda site object and the mechanism to switch between different variants is already implemented. Code to realize my lambda dynamics variant needs to be moved to GROMACS and robbed of C++11 features. Data structures to store and communicate the energy terms (also usable on GPUs) are implemented.

We are very interested in having a general and solid implementation that suites the needs of all of us.

So please let us know if there are additional requirements/things that I missed.

#### #3 Updated by Michael Shirts over 5 years ago

Thomas Ullmann wrote:

OK, I'm going to have to read the paper over more before I comment too much!

For example, the MD integrator should be a function object without any internal state variables which just accepts input coordinates and forces and returns the updated coordinates without knowing about the detailed nature of the coordinates (I think this has been proposed also by others on the developers list for MD integrators in general, If I remember correctly by David van der Spoel).

Correct, though this has always been viewed as a global approach, rather than a local (per site) one. It's not necessarily clear how the different integrators play with each other -- if one set of coordinates is updating with MC, and another with MD, then there may not be theory for how they interact. Clearly, we would want to be able to very easily have rules for changing from one to another globally.

- It should be possible to let different types of force field terms interconvert/switch separately, e.g., Coulomb and Lennard-Jones interactions.

Agreed, see Levi Naden paper. http://pubs.acs.org/doi/abs/10.1021/ct4009188. Similar to Buelens and Grubmüller, but more flexible in a number of ways; it plays very easily with new architectures, since all of the lambda work is done outside of the inner loops.

- it should be possible to consider whole sets of different transformations/reactions

with a user defined free energy calculation method and target statistical error, that is

ideally without the need for postprocessing the simulation results

Why the "no postprocessing" goal?

Consequently, there would be a lambda variable on top of the lambda variables of the sites that controls the transformation. For clarity, I will distinguish this lambda from the lambda values of the site forms by naming it "w" for weight.

In the simplest case there your free energy calculation considers the free energy difference between two forms, 0 and 1, of a single site, w would be given by

w = lambda_0 = 1 -lambda_1. In the general case, there is no such one-to-one correspondence.

Weight isn't a great name, since for many types of dynamics, there is a weight, which is a probability, associated with each vector of lambda.

#### #4 Updated by Thomas Ullmann over 5 years ago

OK, I'm going to have to read the paper over more before I comment too much!

This variant of lambda dynamics is not yet published. Serena's paper dealt with an earlier version.

Correct, though this has always been viewed as a global approach, rather than a local (per site) one. It's not necessarily clear how the different integrators play with each other -- if one set of coordinates is updating with MC, and another with MD, then there may not be theory for how they interact. Clearly, we would want to be able to very easily have rules for changing from one to another globally.

OK. The modularity of the object will, I think, also make it easy to implement the global setting.

Why the "no postprocessing" goal?

I did not mean to discard, e.g., your method of finding the optimal path for the transformation, but to integrate it into the program flow, instead of doing it in a separate step. This is, because we have large-scale computations in mind, where many free energy differences should ideally be computed without user intervention or intermediate processing steps. The goals are user-friendliness and tractability (from the users point of view) of large scale-computations on Excascale computers.

Weight isn't a great name, since for many types of dynamics, there is a weight, which is a probability, associated with each vector of lambda.

OK, maybe there is a better variant to make the distinction clear.

Agreed, see Levi Naden paper. http://pubs.acs.org/doi/abs/10.1021/ct4009188. Similar to Buelens and Grubmüller, but more flexible in a number of ways; it plays very easily with new architectures, since all of the lambda work is done outside of the inner loops.

I knew the paper already and I read it again thoroughly. I think it should play very well with my lambda-dynamics variant. I'm already quite curious to see how it will work for variable sites in a protein environment.

#### #5 Updated by Michael Shirts over 5 years ago

Internally, just having a large set of lambda arrays should be fine. The sophistication is in the setting up which lambda goes to which variable, which is set up by site ID's, form ID's etc. For MC, then we have a (static) list of values for each lambda.

it COULD be set up in a more object oriented way so that lambdas for each site are hidden, but then we lose some ability to quickly address and combine them all, which is likely to lead to some real code slowdown in the future -- I'd like the complicated object oriented stuff to be in the definitions of the indices.

#### #6 Updated by Michael Shirts over 5 years ago

This variant of lambda dynamics is not yet published. Serena's paper dealt with an earlier version.

OK, I don't think we need to talk about lambda dynamics in this thread - I think it can be moved to the other thread, and we can to some extent abstract that.

OK. The modularity of the object will, I think, also make it easy to implement the global setting.

I think this is where posting the documentation would help provide useful feedback from a wider group.

Why the "no postprocessing" goal?

I did not mean to discard, e.g., your method of finding the optimal path for the transformation, but to integrate it into the program flow, instead of doing it in a separate step. This is, because we have large-scale computations in mind, where many free energy differences should ideally be computed without user intervention or intermediate processing steps. The goals are user-friendliness and tractability (from the users point of view) of large scale-computations on Excascale computers.

Ah! OK. I think the key to this is having tools that easily play well with each other. No code can do everything -- instead, the best way is make it easier for things to commuicate with each other. Trying to make one monolithic code that does everything tends to fail. Rather than no postprocessing, I think the better goal is "all postprocessing done transparent handoffs to compatible tools". Getting GROMACS API will be the key to making this very fast. For now, I think python processing for input/output is fine for things that only need to be done every hour or so -- the key is that if it needs to be done every second, then it needs to be internal. So lambda changes -> must be internal now.

I knew the paper already and I read it again thoroughly. I think it should play very well with my lambda-dynamics variant. I'm already quite curious to see how it will work for variable sites in a protein environment.

Great! Levi is watching this thread and can provide feedback on how to do this - he has some other tasks for now. Once the tabulated potentials are in for Verlet, I think this can be built on top very easily.

#### #7 Updated by Gerrit Code Review Bot over 4 years ago

Gerrit received a related patchset '1' for Issue #1652.

Uploader: Thomas Ullmann (thomas_ullmann@web.de)

Change-Id: I68814357b5959ea691bb508fd13aeb8eba4f1fff

Gerrit URL: https://gerrit.gromacs.org/5204

#### #8 Updated by Gerrit Code Review Bot over 4 years ago

Gerrit received a related patchset '1' for Issue #1652.

Uploader: Thomas Ullmann (thomas_ullmann@web.de)

Change-Id: Ic80a3e9de7c54440a2734a6d798479c452383e5e

Gerrit URL: https://gerrit.gromacs.org/5250

#### #9 Updated by Mark Abraham over 3 years ago

**Target version**deleted ()*5.x*

#### #10 Updated by Gerrit Code Review Bot over 2 years ago

Gerrit received a related patchset '5' for Issue #1652.

Uploader: Thomas Ullmann (thomas_ullmann@web.de)

Change-Id: gromacs~master~Ic80a3e9de7c54440a2734a6d798479c452383e5e

Gerrit URL: https://gerrit.gromacs.org/5250

#### #11 Updated by Gerrit Code Review Bot almost 2 years ago

Gerrit received a related patchset '1' for Issue #1652.

Uploader: Thomas Ullmann (thomas_ullmann@web.de)

Change-Id: gromacs~master~Ib042563a6df9b275a6e349a2b5c7434585ac32ac

Gerrit URL: https://gerrit.gromacs.org/7978