JuliaMolSim / Molly.jl

Molecular simulation in Julia
Other
371 stars 51 forks source link

Add anisotropic MC barostat #136

Closed bondrewd closed 1 year ago

bondrewd commented 1 year ago

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

The new interface for using this barostat in 3D and 2D systems is:

MonteCarloAnisotropicBarostat(
    1.0u"bar", # Pressure in the X-axis
    1.0u"bar", # Pressure in the Y-axis
    1.0u"bar", # Pressure in the Z-axis
    temperature,
    boundary
)

and

MonteCarloAnisotropicBarostat(
    1.0u"bar*nm", # Pressure in the X-axis
    1.0u"bar*nm", # Pressure in the Y-axis
    temperature,
    boundary
)

respectively.

If you set any pressure value as 'nothing', the barostat is uncoupled in that axis. For example, for coupling only the Y-axis, you could do the following:

MonteCarloAnisotropicBarostat(
    nothing, # Keep the box in the X-axis fixed
    1.0u"bar", # Scale box in the Y-axis
    nothing, # Keep the box in the Z-axis fixed
    temperature,
    boundary
)

After discussing and merging this PR, I will start working on the constant surface tension MC barostat for membrane simulations.

codecov[bot] commented 1 year ago

Codecov Report

Patch coverage: 92.98% and project coverage change: +0.16 :tada:

Comparison is base (38dec32) 72.90% compared to head (e0425c8) 73.06%.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #136 +/- ## ========================================== + Coverage 72.90% 73.06% +0.16% ========================================== Files 34 34 Lines 4831 4886 +55 ========================================== + Hits 3522 3570 +48 - Misses 1309 1316 +7 ``` | [Impacted Files](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/136?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/136?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2NvdXBsaW5nLmps) | `87.59% <92.45%> (-0.57%)` | :arrow_down: | | [src/spatial.jl](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/136?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL3NwYXRpYWwuamw=) | `76.36% <100.00%> (+0.14%)` | :arrow_up: |

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

jgreener64 commented 1 year ago

Nice one! Something to consider is constructing it with a SVector instead, e.g.

MonteCarloAnisotropicBarostat(
    SVector(1.0u"bar", 1.0u"bar", 1.0u"bar"),
    temperature,
    boundary
)

This gives a little more consistency between the 3D/2D calls and may reduce code reuse internally by allowing one constructor to be defined.

bondrewd commented 1 year ago

Thank you for the comment! I have refactored the constructor, using a SVector for the pressure. It helps to simplify the interface since the dimension of the pressure vector becomes available.

I also added some tests following the current examples. They are designed to check for regressions. In the future, implementing more rigorous tests might be of interest.

jgreener64 commented 1 year ago

The implementation seems good from a quick look.

As you say there is the question of how to test this since we don't currently allow logging pressure in one axis. It might be nice to test with a pressure other than SVector(1.0u"bar", 1.0u"bar", 1.0u"bar") though. Maybe you could test with SVector(1.5u"bar", 0.5u"bar", 1.0u"bar") or similar, and update the ranges for testing as appropriate. Perhaps the box size in each dimension is available to test and shows a difference in this scheme too?

bondrewd commented 1 year ago

Thank you for the suggestion! I added more cases for testing the pressure combination you suggested and different couplings: only X-axis, only Y-axis, only Z-axis, and uncoupled.

jgreener64 commented 1 year ago

Great, thanks.