Started Muller-Brown potential

[ejmeitz]

Hello,

Just started working on the Muller-Brown (MB) potential as a first project so this is definitely a draft. I had some questions about how this potential should be incorporated into Molly. MB is not like any of the other potentials in Molly at the moment, it is only 2D and is not a pairwise interaction (just a surface). Assuming I can figure out how to integrate MB into Molly I will add test coverage in a future commit.

Questions:

1. Other interactions are typed as SpecificInteraction or PairwiseInteraction, MB potential does not fit into either of these. Would I need to create something new (e.g. PotentialEnergySurface)
2. Because it is only 2D how do I make sure this potential is operable within the framework of Molly that is mostly built around pairwise interactions or bonded interactions (e.g. 2 atoms get passed to the potential() function). I have also not implemented a force() function yet. I do not really see how to make the interface used by other potentials amenable to MB.
3. I'm pretty new to Julia so I have a lot of uncertainty around the typing and in general how you go about a large code base without OOP. My brain is stuck thinking in abstract classes and interfaces so any pointers are appreciated.

Thanks, Ethan

[codecov[bot]]

Codecov Report

Base: 83.57% // Head: 83.63% // Increases project coverage by +0.06% :tada:

Coverage data is based on head (45c43a2) compared to base (b172d0d). Patch coverage: 90.90% of modified lines in pull request are covered.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #105 +/- ## ========================================== + Coverage 83.57% 83.63% +0.06% ========================================== Files 32 33 +1 Lines 3903 3936 +33 ========================================== + Hits 3262 3292 +30 - Misses 641 644 +3 ``` | [Impacted Files](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/105?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim) | Coverage Δ | | |---|---|---| | [src/Molly.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/105?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL01vbGx5Lmps) | `100.00% <ø> (ø)` | | | [src/interactions/muller\_brown.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/105?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2ludGVyYWN0aW9ucy9tdWxsZXJfYnJvd24uamw=) | `90.90% <90.90%> (ø)` | | Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim)

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

[JaydevSR]

Hi @ejmeitz, I think you are looking to implement something referred to as a General Interaction in Molly's framework. One example is: https://github.com/JuliaMolSim/Molly.jl/blob/9ced8f7261807467f928b4ff207f0fac96f4c525/src/interactions/implicit_solvent.jl#L338

I think a potential surface can be implemented similarly. As for the typing, the basic idea is to use variable type parameters in the struct definitions and then in the respective constructor we figure out if the types of the inputs to the constructor make sense and if they do we initialize the struct with those parameters and inputs. If this is done right everything works together nicely. You can see this in all the types (structs) defined in Molly. Take a look at the file src/types.jl for numerous examples.

I think @jgreener64 will be able to clarify further if I am wrong.

[jgreener64]

Thanks for giving this a go @ejmeitz.

@JaydevSR is correct, this can be implemented as a general interaction (https://juliamolsim.github.io/Molly.jl/dev/docs/#General-interactions). In the first instance you can write a function that takes in a 2D coordinate and gives the force, and also one for the energy. These functions can then be broadcasted over the system coordinates inside the forces and potential_energy functions. Broadcasting works here since the force/energy only depends on the particle being considered.

[ejmeitz]

@JaydevSR @jgreener64 thank you for the tips! I think I have the code conforming to the framework correctly now. I'm working on tests, but still had a few questions about the code I wrote.

1. Should I enforce that an input has certain units? e.g. the A vector should have units of energy, x0 units of length etc. Or is that on the user to keep consistent?
2. Should I enforce that the user only inputs 2D data into this system? The broadcasting in the forces and potential_energy is agnostic of the dimension of each sub-array currently, but should only work for 2D.
3. Formats for output of forces() and potential_energy()
   • Should PE be a single number for whole system or a single number for each atom
   • Force should have info for every particle but does the type of the output matter as long as its some kind of nested array?

I can probably answer 3 myself, but I am still setting up my dev environment for Molly to do local testing.

[jgreener64]

Should I enforce that an input has certain units?

Leave it to the user to ensure consistency, this is usually fine since arithmetic with inconsistent units will throw an error.

Should I enforce that the user only inputs 2D data into this system?

If you replace coord with coord::SVector{2} in the arguments of potential then an error will be thrown in non-2D cases.

Should PE be a single number for whole system or a single number for each atom

For the whole system.

Force should have info for every particle but does the type of the output matter as long as its some kind of nested array?

Should be SVector{2, T} for each atom (T is determined by the force unit).

[ejmeitz]

@jgreener64 pretty sure the code is all correct now, but I can't figure out how to get the codecov tests to pass. I added some tests and they cause the ubuntu check to fail. For now they're just commented out at the bottom of interactions.jl in my last commit.

I also tried plotting the trajectories as a sanity check for myself to see if they match this Wolfram notebook. In my simulation the particle just goes to the bottom of the energy potential (see below) instead of oscillating around the minimum which was a little un-expected. The energy contours and forces looked alright so I'm not sure why I did not get the correct dynamics (used VV integrator for a few million steps).

Thanks for your patience on this first pull-request! Once I figure out this last piece I plan on starting the Gay-Berne and NEQUIP potentials or the NoseHoover thermostat.

[jgreener64]

You can run the tests locally or check the logs (https://github.com/JuliaMolSim/Molly.jl/actions/runs/3833098088/jobs/6524174032) to see the issue. In this case testing that force(inter, local_min2) is approximately equal to 0.0u"kJ * mol^-1 * nm^-1" fails since the first is an SVector{2} and the second is a single value.

In my simulation the particle just goes to the bottom of the energy potential (see below) instead of oscillating around the minimum which was a little un-expected.

This will depend on the simulator you use, I don't know which one the Wolfram notebook used. You would see more oscillation with a Berendsen thermostat for instance.

The code looks broadly okay, I'll add some review comments once tests pass.

[ejmeitz]

Couldn't get tests for "potential" to work as that function is not exposed to the test set (or maybe I'm being dumb). Can also write an example for docs if necessary.

[ejmeitz]

@jgreener64 Renaming the force function to force_muller_brown breaks something in the simulation loop. Molly is trying to call force() but finds nothing. If I leave it as force things seem to run fine. It does not seem to care what I call the potential energy function. All tests pass if I do not rename but I want to make sure I understand how everything is interacting here. Like you said I thought forces() and potential_energy() would be called not force() and potential().

Thanks, Ethan

[jgreener64]

If you rename force to force_muller_brown then simulation should not be calling force at all, so look at the line number in the error stacktrace and see where that is. Maybe you forgot to change where the function is broadcasted in forces? Potential energy is not called during simulation unless you use an energy logger.

[ejmeitz]

Any idea why the minimization test fails on my laptop but passed in the GitHub checks? I get an error on line 62: @test all(x -> isapprox(x, 0.4u"nm"; atol=1e-3u"nm"), dists_flat) which is kinda weird cause I never touch that code and it passes on github. I think the tests I added should pass but am unable to get past the minimization test on my laptop for now....will try again on desktop later.

[jgreener64]

Does it still fail on repeated runs? There is some stochasticity in the tests.

It could also be something with GPU setup since that is a GPU test, what OS and GPU are you on?

You could also print out the distance values and see what is going on.

[ejmeitz]

Yeah I ran is like 4 times last night and printed them they are all within the expected tolerance. Here's what I get re-running this morning: [0.4000429540445738 nm, 0.4004080941145199 nm, 0.40004454772593906 nm] and it still fails even though the expected val is 0.4 and the atol is 0.001

Using intel i7-10510U, integrated graphics and there's an NVIDIA MX250 also which is what the tests run on I'd imagine.

[jgreener64]

That's strange, does this return true or false? What about without the unit annotations?

using CUDA, Unitful
dists_flat = CuArray([0.4000429540445738, 0.4004080941145199, 0.40004454772593906]u"nm")
all(x -> isapprox(x, 0.4u"nm"; atol=1e-3u"nm"), dists_flat)

[ejmeitz]

It returns true with and without units. Super weird, tests passed on GitHub though.

This is full error, don't think I'm missing anything:

Energy minimization: Test Failed at C:\Users\ejmei\Repositories\Molly.jl\test\minimization.jl:62
  Expression: all((x->begin
            isapprox(x, 0.4 * u"nm"; atol = 0.001 * u"nm")
        end), dists_flat)
 [4] top-level scope
   @ C:\Users\ejmei\Repositories\Molly.jl\test\minimization.jl:2
Test Summary:       | Pass  Fail  Total     Time
Energy minimization |    5     1      6  1m43.9s
ERROR: LoadError: Some tests did not pass: 5 passed, 1 failed, 0 errored, 0 broken.
in expression starting at C:\Users\ejmei\Repositories\Molly.jl\test\minimization.jl:1
in expression starting at C:\Users\ejmei\Repositories\Molly.jl\test\runtests.jl:61
ERROR: Package Molly errored during testing

[jgreener64]

Let's just call it weirdness for now.

Thanks for doing this PR!