Constant pH module
===

###### tags: `constant pH` `Meeting`

:::info
- **Location:** [zoom-link](https://kth-se.zoom.us/j/62870145132)
- **Date:** 2022-05-04 11:00
:::

## Notes on the current code

- For multidimensional matrices like `constraintMultipliers` use `MultiDimArray` from `multidimarray.h` like `MultiDimArray<std::vector<float>, dynamicExtents3D>> constraintMultipliers(n_rows, n_columns)`

- Shouldn't the constant-pH modules be implemented to be used with all integrators (not just md?) Especially energy minimsation is something that I believe we want.
- This will require that all things in `md.cpp` will have to move down to `do_force` or up to `runner.cpp`

- -

## Design for charge setter

- `ChargeSetterManager` lives in MdAtoms
- modules need access to `ChargeSetterManager` before `do_force` loop and after `LocalAtomSets` are established -> add `setupNotification` in `runner.cpp`
- `ChargeSetterManager` holds a const view on `chargeB`, a copy of `chargeA` and a view on `chargeA` to modify
- `ChargeSetterManager` makes sure that there is no overlap between atom indices trying to set charges
- `ChargeSetterManager` obtains localatomsets
- `ChargeSetterManager` should be instantiated after localatomsetmanager

digraph chargesetter {
    pad= 2.0
    splines= ortho
    nodesep= 0.60
    ranksep= 0.75
    fontname= "Sans-Serif"
    fontsize= 13
    fontcolor= "#2D3436"

    node [shape = box
        // style = rounded
        fixedsize = true
        width = 1.4
        height = 1.4
        labelloc = c
        fontname = "Sans-Serif"
        fontsize = 13
        fontcolor = "#2D4436"]
        
    atom_indices_1 [color = salmon, label="atom\nindices\n1"]
    atom_indices_2 [color = lightblue, label="atom\nindices\n2"]
    chargeA [label = "chargeA\ntopology"]
    chargeB [label = "chargeB\ntopology"]
    chargesettermanager [label = "charge\nsetter\nmanager"]
    lambda_1 [color = salmon, label = "lambda 1"]
    lambda_2 [color = lightblue, label = "lambda 2"]
    chargeAsim [label = "chargeA\nsimulation"]
    chargesetter_1 [color = salmon, label = "charge\nsetter\n1"]
    chargesetter_2 [color = lightblue, label = "charge\nsetter\n2"]
    atom_indices_1 -> chargesettermanager [color = salmon]
    atom_indices_2 -> chargesettermanager [color = lightblue]
    chargeA -> chargesettermanager
    chargeB -> chargesettermanager
    chargesettermanager -> chargeAsim
    lambda_1 -> chargesetter_1 [color = salmon]
    lambda_2 -> chargesetter_2 [color = lightblue]
    
    chargesettermanager -> chargesetter_1 [color = salmon]
    chargesettermanager -> chargesetter_2 [color = lightblue]
}

### Questions

- do we have to set the charges for all global atoms on all nodes? (I would guess so)

- for masses we use a separate `md.tmass` and seperate that out from massA and massB etc, but for charges we use `chargeA` both as the *default* charge and as basis for interpolation with lambda. Then the kernels do the interpolation in the interaction functions, but is there really any gain in this? Wouldn't it be much simpler to also have a `charge` array that is interpolated in all FE steps?

## Design for potential calculator

- add `vector<bool>` is in set to localatomset if requested
### Current state

- currently `fr` is accessbile in kernel
- in listed-forces `pairs.cpp` access `fr->electrostaticPotential` which is a view on an array that is owned by the constant pH module.
- requires also a `vector<bool>` to speed up query whether an atom is a lambda atom

### Questions
- How do you treat the long-range contributions to the potential?

## Integration trigger

- add callback when do_force is done, so that we know all potential calculations are ready. Then lambda will be updated and charges set.

## Data flow and residency

digraph structs {
    pad= 2.0
    splines= ortho
    nodesep= 0.60
    ranksep= 0.75
    fontname= "Sans-Serif"
    fontsize= 13
    fontcolor= "#2D3436"

    node [shape = box
        // style = rounded
        fixedsize = true
        width = 1.4
        height = 1.4
        labelloc = c
        fontname = "Sans-Serif"
        fontsize = 13
        fontcolor = "#2D4436"]
    subgraph clustercpH {
        
        style=filled
        penwidth=0
        fillcolor = "#FFEEEE"
        d_lambda [label = "lambda'"]
        lambda
        atom_indices [label = "c-pH\n atom indices"]
        lambda_update [shape = circle, label = "lambda\nupdate"]
        potential_retreiver [label = "potential\nretreiver"]

    }
	ckpt [shape = circle, label="checkpointing"]
	charge_update [label="charge\nupdater"]
    potential_calculator [label = "potential\ncalculator"]
    force [shape = circle, label="force\ncalculation"]
        lambda->ckpt
        d_lambda -> ckpt
        ckpt -> lambda
        ckpt -> d_lambda
        lambda -> charge_update
        charge_update -> force
        atom_indices -> potential_calculator
        potential_calculator -> potential_retreiver
        potential_retreiver -> lambda_update
        force -> potential_calculator
        lambda_update -> lambda
        lambda_update -> d_lambda
}

## Questions
How do you update the charges now? I would do it in do_force, but there we would have to remove the const qualifier from mdatoms a long way back. The alternative would be to hand down an object that can modify mdatoms, but this is also not so nice, because you could be tricked into thinking that mdatoms are const. This charge-modifier would live in MdAtoms then.

void do_force(FILE*                               fplog,
              ...
              tensor                              vir_force,
              const t_mdatoms*                    mdatoms,
              const MdAtomsChargeModifier &       chargeModifier,
              gmx_enerdata_t*                     enerd,
              ...

https://gitlab.com/gromacs/gromacs/blob/21d738a681a62f18fc795e57dabddb4f9a0708e2/src/gromacs/mdlib/sim_util.cpp#L1316

---

What should we focus on now - is the multi-topology implementation within reach? Should we focus on finishing the constant pH module with the hack of modifying the charges above?

---

What about all other types of interaction terms? Change of atom type, etc..

Core idea MdModules - callbacks vs interface classes
----------------------------------------------------

Desing principles MdModules

1. only modules can modify their data
- keep code in one place
- no surprises
- easier to test - no need to call routines from all over the codebase to test a module
- no need to understand the internal workings of other parts of the code
- enforces writing clean interfaces
-