implement MC membrane barostat

[bondrewd]

This PR adds an implementation for the Monte Carlo membrane barostat following OpenMM's implementation.

This barostat is only available for 3D systems, and the API is as follows:

MonteCarloMembraneBarostat(
    1.0u"bar",   # Pressure in the system
    1.0u"bar*nm" # Surface tension
    temperature,
    boundary;
    xy_isotropy=false,        # If true, scale the XY-axes isotropically
    constant_volume=false,    # If true, keep the volume constant by scaling the Z-axis accordingly
    z_axis_fixed=false        # If true, keep the Z-axis fixed
)
and

In this barostat, constant_volume and z_axis_fixed cannot be true simultaneously.

[codecov[bot]]

Codecov Report

Patch coverage: 94.02% and project coverage change: -8.69 :warning:

Comparison is base (d7a0d90) 73.04% compared to head (ce3fb4d) 64.36%.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #139 +/- ## ========================================== - Coverage 73.04% 64.36% -8.69% ========================================== Files 34 34 Lines 4886 4944 +58 ========================================== - Hits 3569 3182 -387 - Misses 1317 1762 +445 ``` | [Impacted Files](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/139?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim) | Coverage Δ | | |---|---|---| | [src/coupling.jl](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/139?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2NvdXBsaW5nLmps) | `87.75% <94.02%> (+0.83%)` | :arrow_up: | ... and [13 files with indirect coverage changes](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/139/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim)

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

[bondrewd]

Sure! I will open an issue to discuss the example and related contributions.

[jgreener64]

There seem to be occasional test failures with this barostat for the constant_volume=true case. See for example https://github.com/JuliaMolSim/Molly.jl/actions/runs/5466737179/jobs/9952101066.

Any chance you could take a look?

[bondrewd]

Sure! Today I tested the MC membrane barostat and found that for some seeds, the 25-atoms-system adopts conformations in which one axis is favored.

I tried increasing the number of atoms, e.g., 100, and found the same behavior:

Overall, it seems that an LJ gas might not be the most suitable system for testing the MC membrane barostat. The reason is that the surface tension term \gamma * dA favors the scaling in the XY plane (for a positive surface tension). For membrane systems, there exists a symmetry in the X- and Y-axes that balance the scaling. For a gas that is fully symmetric (point masses with radial LJ potential), there is no counterbalance for the surface tension.

One possible solution is to decrease the scaling factor, the scaling increment, and the volume frac. This effectively weakens the coupling, increasing the stability of the system. However, I am not quite sure if this is the best approach since it feels like a band-aid solution.

Another idea is to momentarily comment out the test until the MARTINI PR is merged. Then we could, for example, use a small and equilibrated membrane system composed of, let's say, 128 lipids, and @test some properties such as the pressure, the volume, the area per lipid, and the membrane thickness. Since a coarse-grained run is relatively cheap to compute, it might not have a huge impact on the CI execution time. @jgreener64, what do think about this approach?

The other issue that needs to be addressed is the implementation of correctness tests since what I have described above only covers regressions. One possible way is to sample with the barostat and test the fluctuations in volume and pressure. This might be a bit off-topic, but testing MD code is not so trivial and we might want to study how other more mature software are doing it.

[jgreener64]

Thanks for looking into it! Yes in the short term I will just remove the tests in if barostat.constant_volume and long term we can use MARTINI to test pressure more thoroughly.

I agree that correctness testing is difficult, especially when stochasticity is involved. I have consistently run into this when implementing algorithms for the package. In general it would be useful to study other software and see how they test correctness, or at least try to match them exactly for more algorithms (currently we match OpenMM exactly for a short NVE protein simulation). Consistent random numbers across languages is one issue there.