---
# System prepended metadata

title: Constant pH implementation in GROMACS

---

---
type: slide
slideOptions:
  theme: white
---

# Constant pH implementation in GROMACS
## Additional functionality required from MDModules 

---

## Brief review of algorithms

- Protonation is treated as independent $\lambda$-coordinates, propagated alongside with cartesian coordinates. Point charges are defined by $\lambda$-values: $q_i = (1-\lambda)q_i^{\text{prot}}+\lambda q_i^{\text{deprot}}$
- Forces on $\lambda$-coordinates is computed as $\sum\Delta q_i \phi(x_i)$, where sum is taken over all atoms that constitute $\lambda$-group, and $\phi(x_i)$ - is the potential at atom $i$ with all charges interpolated
- Additional biasing forces on $\lambda$-coordinates are required. They depend only on $\lambda$-values

---

## Brief review of algorithms

- Box neutrality is maintained by charge constraints. Multiple ($N_{\text{buf}}$) weakly interacting titratable buffers are introduced into the system and total charge ($\sum q_{\lambda_i} + N_{\text{buf}} \cdot q_{\text{buf}}$) is constrained, by constraining the corresponding linear combination of $\lambda$-values
- Multisite residues (e.g. Histidine) are represented by several $\lambda$-coordinates (each coordinate correspond to a physical states). The sum of those coordinates is constraied to 1.

---

## Implementation details

We plan to implement constant pH using MDmodules

```sequence
    participant MainLoop
    LambdaDynamics->>MainLoop: Set charges
    MainLoop->>MainLoop: Compute forces and potentials
    MainLoop->>LambdaDynamics: Send potentials
    LambdaDynamics->>LambdaDynamics: Update and constrain lambdas
```

---

## Structure of IMDmodule

```cpp=
class IMDModule
{
public:
    virtual ~IMDModule() {}

    //! Returns an interface for handling mdp input (and tpr I/O).
    virtual IMdpOptionProvider* mdpOptionProvider() = 0;
    //! Returns an interface for handling output files during simulation.
    virtual IMDOutputProvider* outputProvider() = 0;
    //! Initializes force providers from this module.
    virtual void initForceProviders(ForceProviders* forceProviders) = 0;
    //! Subscribe to simulation setup notifications
    virtual void subscribeToSimulationSetupNotifications(MDModulesNotifiers* notifiers) = 0;
    //! Subscribe to pre processing notifications
    virtual void subscribeToPreProcessingNotifications(MDModulesNotifiers* notifiers) = 0;
};
```

---

## Additional callback

- Constant pH actually doesn't provide any forces on atoms, thus it is better to not perform $\lambda$-dynamics loop in `ForceProvider`
- Idea is to add additional callback `MDLoopNotification` in parallel with `ForceProvider`
- Question: Do we want to move parts `QMMM` and `DensityFitting` force providers, so only forces are computed in module `ForceProvider` and helper jobs are done in `MDLoopNotification`?


---

## Structure of IMDmodule option provider

```cpp=
class IMdpOptionProvider
{
public:
    virtual void initMdpTransform(IKeyValueTreeTransformRules* transform) = 0;
    virtual void initMdpOptions(IOptionsContainerWithSections* options) = 0;
    virtual void buildMdpOutput(KeyValueTreeObjectBuilder* builder) const = 0;
};
```

---

## Structure of MDModules

```cpp=
class MDModules
{
public:
    MDModules();
    ~MDModules();

    void initMdpTransform(IKeyValueTreeTransformRules* rules);
    void buildMdpOutput(KeyValueTreeObjectBuilder* builder);
    void assignOptionsToModules(const KeyValueTreeObject& params, IKeyValueTreeErrorHandler* errorHandler);
    void adjustInputrecBasedOnModules(t_inputrec* ir);
    IMDOutputProvider* outputProvider();
    ForceProviders* initForceProviders();
    void subscribeToSimulationSetupNotifications();
    void subscribeToPreProcessingNotifications();
    void add(std::shared_ptr<IMDModule> module);
    const MDModulesNotifiers& notifiers();
};
```

---

## Reading options for MDModules

1. Initialization of `mdp` options (all left side names have to be known)

2. Assigning `mdp` options to modules

```cpp
mdModules->initMdpTransform(transform.rules());
```
```cpp!
mdModules->assignOptionsToModules(*ir->params, &errorHandler);
```

---

## Input organisation

Lambda Dynamics needs three types of input:

- General input (lambda dynamics is on, lambda thermostat parameters, lambda particles properties, etc.)
- Description of $\lambda$-group types (how many states does the group have, which atom types build the group, what are internal FF-dependent parameters of each state (e.g. pKa, dvdl-coefficients))
- Mapping of index groups to $\lambda$-group types

---

## Mapping of index groups to group types

Atomic indices are stored in index groups and are mapped to group types using `.mdp` file (similar to pulling groups)

```
    lambda-dynamics-atom-set1-name                = ASP
    lambda-dynamics-atom-set1-index-group-name    = LAMBDA1
    lambda-dynamics-atom-set1-barrier             = 5.0 
    lambda-dynamics-atom-set1-initial-lambda      = 1.0
```

---

## Additional routines for reading dependent `mdp` options

```cpp
mdModules->initMdpDependentTransform(
    IKeyValueTreeTransformRules* rules, 
    const size_t nGroups);

mdModules->assignOptionsToModules(
    const KeyValueTreeObject& params, 
    IKeyValueTreeErrorHandler* errorHandler,
    const size_t nGroups);
```

You can comment on the discussion in MR [!3135](https://gitlab.com/gromacs/gromacs/-/merge_requests/3135)