Monte-Carlo for molecules

[JaydevSR]

This PR adds some simple Monte-Carlo functionality i.e., Metropolis Monte-Carlo for molecules. I was initially thinking of writing the simulator in simulators.jl but there may be some monte-carlo specific code (like the random translation trial moves here) that doesn't fit there so I added a new file.

Some things that need to be done before merge:

• [x] Use of neighbor finder for calculating potential energy at each step?
• [ ] Each step uses copy() which is not very efficient, any suggestions for improvement?
• [x] Example simulations (one example given below)
• [x] Docs and Tests

I have run a test simulation given below:

n_atoms = 100
atom_mass = 10.0u"u"
atoms = [Atom(mass=atom_mass, charge=10, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]
boundary = CubicBoundary(2.0u"nm", 2.0u"nm", 2.0u"nm")
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm")
pairwise_inters = (Coulomb(),)

temp_mult = [1598.0u"K", 1198.0u"K", 798.0u"K", 398.0u"K", 198.0u"K", 98.0u"K", 8.0u"K"]
sys = System(
      atoms=atoms,
      pairwise_inters=pairwise_inters,
      coords=coords,
      boundary=boundary,
      loggers=(coords=CoordinateLogger(n_atoms),),
)

trial_args = Dict(:shift_scaling => 0.1, :length_units => u"nm")
for t in temp_mult
    sim = MetropolisMonteCarlo(; 
            temperature=t,
            trial_moves=random_normal_translation!,
            trial_args=trial_args
        )

    simulate!(sys, sim, 20000)
end

This is what the simulation looks like:

https://user-images.githubusercontent.com/59474411/199914858-3d8d077e-9701-4ad5-abcd-48c06018acc4.mp4

As expected the atoms are arranging into a lattice (sort of). Another simulation that was suggested by @jgreener64 was to sample the torsional angle of ethane. I am not quite sure how to set that up so some help will be appreciated.

[codecov[bot]]

Codecov Report

Base: 83.87% // Head: 83.74% // Decreases project coverage by -0.13% :warning:

Coverage data is based on head (93d16bb) compared to base (47f3e1f). Patch coverage: 84.78% of modified lines in pull request are covered.

:exclamation: Current head 93d16bb differs from pull request most recent head ab4b6e6. Consider uploading reports for the commit ab4b6e6 to get more accurate results

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #104 +/- ## ========================================== - Coverage 83.87% 83.74% -0.14% ========================================== Files 32 33 +1 Lines 3796 3844 +48 ========================================== + Hits 3184 3219 +35 - Misses 612 625 +13 ``` | [Impacted Files](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/104?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/104/diff?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/montecarlo.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL21vbnRlY2FybG8uamw=) | `81.08% <81.08%> (ø)` | | | [src/loggers.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2xvZ2dlcnMuamw=) | `95.65% <100.00%> (+0.19%)` | :arrow_up: | | [src/zygote.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL3p5Z290ZS5qbA==) | `79.61% <0.00%> (-0.94%)` | :arrow_down: | | [src/interactions/implicit\_solvent.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2ludGVyYWN0aW9ucy9pbXBsaWNpdF9zb2x2ZW50Lmps) | `90.00% <0.00%> (-0.07%)` | :arrow_down: | 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.

[jgreener64]

Looks like a great start. If the lattice formation is clearer in 2D you could do that. I think the copy() is okay for the moment, profiling would help point to the slow steps though.

Use of neighbor finder for calculating potential energy at each step?

I think call find_neighbors every step and use the neighbours, as the other simulators do. In many cases, for example with large trial moves, not using a neighbour list would be more efficient since neighbour calculation would need to be run every step. It can be left to the user though since the default is NoNeighborFinder.

sample the torsional angle of ethane

I can look at this myself after the PR is merged if it's easier.

A couple of other points:

• The potential_energy only needs to be calculated once per step, with E_old being updated at the end of each step.
• It would also be nice to have a logger similar to the ReplicaExchangeLogger that logs energies and acceptance/failure.
• For the random translation functions I would make the norm of the displacement match the stated uniform/Gaussian distribution, rather than doing it separately by dimension. This can be done by generating x/y/z displacements with randn, dividing by the norm to get a unit vector, and multiplying by the distance drawn from the appropriate distribution.

[JaydevSR]

It would also be nice to have a logger similar to the ReplicaExchangeLogger that logs energies and acceptance/failure.

So, for the logger I am thinking that instead of logging the energies for which we already have a logger, we can log the value of $\displaystyle\frac{E}{k_B T}$ i.e. the argument of the Boltzmann rates which makes sense for a Monte-Carlo simulation. What do you think? Also, should I just log the number of acceptances and selections or also a bit array of when the acceptances occur?

[JaydevSR]

A 2D simulation which shows lattice formation with 100 atoms in 4nm square boundary:

https://user-images.githubusercontent.com/59474411/200545210-7b8d61c5-e23a-445a-8eef-eda561721dd0.mp4

Similarly for 3D with 4nm boundary and 100 atoms:

https://user-images.githubusercontent.com/59474411/200545464-c5f4fffe-c61d-4933-a8db-a0df5a1383b9.mp4

[jgreener64]

Yes I think logging that value sounds good. I would log the bit array of when the acceptances occur too.

[JaydevSR]

Now the question remains is where the logger would go. If we put it in loggers named tuple in the System type, then we cannot pass in info about whether the acceptance occurred or not. And also, it would recalculate any potential energy that we have already calculated in the metropolis part. In the replica exchange, I added a field to replica system. So, should there be a similar optional thing for the System type, like a montecarlo_logger field?

[jgreener64]

This seems to be a general problem when writing simulator-specific loggers. I think we should change the arguments of the run_loggers! and log_property! function definitions to accept additional kwargs.... In this case we could then pass the acceptance info and potential energy as keyword arguments to run_loggers! and use them in the relevant log_property! function. I think this should mostly be a case of adding kwargs... to the keyword arguments of functions in src/loggers.jl.

[JaydevSR]

I have added all of the things that are needed for the PR ✌️. Please tell me if some more tweaks are needed or if I missed something.

[JaydevSR]

I have added tests for acceptance rate (test for bad cases), avg potential energy (should decrease from a random config at start to end for a repulsive system), and the mean nearest neighbor distance (for a repulsive system it should not be less than the Wigner-Seitz radius and for the given boundary it is not possible for it to be 2 times the Wigner-Seitz radius).

[jgreener64]

Great, thanks for this.