NoseHoover Thermostat

[ejmeitz]

Added NoseHoover thermostat to the simulations file. The code is 95% done and I will add tests soon.

Some questions I have:

• I do not think this could be created as a "coupler" due to the presence of zeta which needs to be updated in the simulation loop, am I missing something?

Unrelated to This Pull Request:

• Is the neighbor list being updated every timestep (looks like it)?
• In the code below why are constraints applied to the copied coordinates? Those coordinates are never used again and not updated by the Integrator.

old_coords = copy(sys.coords)
 sys.coords += sys.velocities .* sim.dt .+ ( (remove_molar.(accels_t) .- zeta.*sys.velocities).* sim.dt ^ 2) ./ 2

 apply_constraints!(sys, old_coords, sim.dt)
 sys.coords = wrap_coords.(sys.coords, (sys.boundary,))

[jgreener64]

Looks good so far. Out of interest are you basing this on a particular implementation in another software? I've seen a few different ways to implement Nosé–Hoover around.

I do not think this could be created as a "coupler" due to the presence of zeta which needs to be updated in the simulation loop

Yes I agree.

Is the neighbor list being updated every timestep

Not in general, no. find_neighbors is called every time step and then it is up to the neighbor list what to do. Most have a property n_steps that determines the number of steps between updates. If it is not time to update then the current neighbors are returned, see the !iszero(step_n % nf.n_steps) && return current_neighbors lines in the neighbors.jl file. Calling find_neighbors every step keeps flexibility since it means that these decisions are handled by the neighbor list implementation.

In the code below why are constraints applied to the copied coordinates?

Applying the constraints requires the old coordinates old_coords and the updated coordinates sys.coords, so both are given to apply_constraints!.

[ejmeitz]

It should just be the same as velocity verlet (it matched your code for VV). The ODE being solved is below but that should be universal. If you're referring to something else let me know I've only seen this implementation. I can check what LAMMPS uses if you're still curious. r_dot = v v_dot = F/m - zeta * v zeta_dot = (1/tau^2)(T/T_des - 1)

[ejmeitz]

Assuming the tests past after this latest commit, I think it should be done.

Do not really know how to write tests for this, temperature will fluctuate around target and the energy should follow a normal distribution but there's not a super rigorous way to check that either of these are true.

[codecov[bot]]

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.13 :tada:

Comparison is base (fee9cac) 83.63% compared to head (7fad457) 83.77%.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #107 +/- ## ========================================== + Coverage 83.63% 83.77% +0.13% ========================================== Files 33 33 Lines 3936 3969 +33 ========================================== + Hits 3292 3325 +33 Misses 644 644 ``` | [Impacted Files](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/107?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim) | Coverage Δ | | |---|---|---| | [src/simulators.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/107?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL3NpbXVsYXRvcnMuamw=) | `95.95% <100.00%> (+0.39%)` | :arrow_up: | | [src/zygote.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/107?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL3p5Z290ZS5qbA==) | `79.53% <0.00%> (ø)` | | | [src/interactions/muller\_brown.jl](https://codecov.io/gh/JuliaMolSim/Molly.jl/pull/107?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL2ludGVyYWN0aW9ucy9tdWxsZXJfYnJvd24uamw=) | `90.90% <0.00%> (ø)` | | 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]

For testing you could run a simulation (e.g. Lennard-Jones) and test that the mean/std of temperature and energy are roughly as expected. Doesn't need to be too rigorous but it will mean that the simulator is run in the tests and prevent obvious regressions. Might be good to see some plots of e.g. temperature over the simulation too.

[ejmeitz]

Triple checked integration, units and system set up...cannot figure out why things blow up. Will come back to later.

[jgreener64]

Don't know if this is the issue but you should make min_dist for place_atoms larger than the LJ σ since the minimum of the potential is at 2^(1/6) * σ.

[ejmeitz]

That did not change much unfortunately. The system setup is also what I think is incorrect because I have no way of knowing if what I initialize is a gas, liquid or solid or something in between.

Do you know of a way to instantiate particles in FCC/BCC instead of placing randomly? Pretty sure I can load an external file, but would be nice to generate through code. I have this on a list of things I was going to add, but have not looked at the external libraries Molly imports to see if this feature exists already.

[jgreener64]

That would be a nice addition. There may be something relevant in https://github.com/mdavezac/Crystals.jl.

You could also try using the protein from the example at https://juliamolsim.github.io/Molly.jl/stable/docs/#Simulating-a-protein by switching the simulator, this should be stable.

[ejmeitz]

Protein system still blows up temperature to infinity, guess structure isn't the problem (or only problem).

Crystals.jl looks interesting, I'll definitely work on that as my next project.

[ehsanirani]

I have no way of knowing if what I initialize is a gas, liquid or solid, or something in between.

One quick comment regarding the initial configuration: I usually start with a soft potential instead of LJ, then minimize the potential energy, so overlaps are naturally removed. Then I replace the soft potential with LJ and run the main simulation.

[ejmeitz]

Yeah I'm just gonna put two particles next to eachother so that they oscillate. I also have code from a class for NoseHoover that should just spit out the same thing. Should be fixed soon.

[ejmeitz]

Fixed bug... still writing tests and trying to convince myself its all OK by calculating heat capacity from the energy distribution.

[ejmeitz]

Ok I added a test, but I was never able to convince myself the heat capacity I get from the energy distribution is correct. I can get a number from the simulation, but its not immediately clear to me what to compare it to. Will come back ad check this once I add the FCC/BCC stuff.

I also couldn't figure out why the temperature plots look different between this and LAMMPS (ignore the different temp the behavior is independent of that). I don't think its objectively wrong it just seemed weird to me: LAMMPS NoseHoover: Molly NoseHoover:

[jgreener64]

If it's proving hard to test it might be worth trying to set up exactly the same system in another software like LAMMPS and checking that the behaviour is the same.

[ejmeitz]

Will update the test later to remove stochasticity. In the mean time, I ran the protein simulation from the docs with NoseHoover. Temperature and energy plots are pasted below in comparison to the Langevin integrator used by the example.

Langevin:

Nose Hoover

[ejmeitz]

Yeah now that you mention it my equation for zeta looks different than anything I find online. This was just the form I was taught in a class I'll look for a source or try and show that it is equivalent to something from online.

The ODE I was taught has zeta_dot = (1/damping^2)(T_inst/T_desired - 1) which looks different than most of the things I see online.

[ejmeitz]

@jgreener64 Here's a source: https://aip.scitation.org/doi/pdf/10.1063/1.449071

I can upload PDF if needed or put in the comments in the code.

[jgreener64]

I'm struggling to link the equations in the paper to the integrator in the code, I'll admit that my understanding of the Hamiltonian formulation of dynamics is somewhat limited though.

Perhaps @noeblassel could have a look and comment on the implementation, or provide a reference to the integration scheme used?

[ejmeitz]

This paper is more just a formulation that writes the ODE in a way that matches how I integrated it. The numerical scheme is exactly the same as Velocity Verlet (I even copied most of the code from there). The only addition is zeta. The ODE's I integrate with that scheme are written in equation 22. The Q is defined in 24, and makes the zeta equation (1/tau^2)*(T_actual/T_des - 1).

[noeblassel]

Hi, looks pretty good, I had a quick peak at the scheme:

    `v_half = sys.velocities .+ ((sim.dt/2).* (remove_molar.(accels_t) .- (zeta.*sys.velocities)))
    old_coords = copy(sys.coords)
    sys.coords += (sim.dt .* v_half)

    KE_half = sum(masses(sys) .* sum.(abs2, v_half)) / 2
    T_half = uconvert(u"K",2 * KE_half / (df * sys.k))
    zeta += (sim.dt / (sim.damping^2)) * ((T_half/sim.temperature) - 1)

    apply_constraints!(sys, old_coords, sim.dt)
    sys.coords = wrap_coords.(sys.coords, (sys.boundary,))

    accels_t_dt = accelerations(sys, neighbors; n_threads=n_threads)

    sys.velocities = (v_half .+ (sim.dt/2)*(remove_molar.(accels_t_dt)))./(1 + (zeta*sim.dt/2))

    sim.remove_CM_motion && remove_CM_motion!(sys)
    apply_coupling!(sys, sim.coupling, sim)`

I agree it ressembles a Velocity-Verlet like scheme for the ODE described in equation (8) -- with constants factored out into the sim.damping parameter, except for the final velocity update which I imagine should instead be sys.velocities =v_half + (sim.dt/2) * (remove_molar.(accels_t) .- (zeta.*sys.velocities))) -- although this is consistent to the first order in sim.dt with the one proposed.

To really get a symplectic scheme, I think you would need to keep track of the $s$ variable (e.g. equation 9 in the reference), and update the auxilliary state variable $(s,\xi)$ according to a Velocity-Verlet scheme for the Nose Hamiltonian.

Here I think xiis only updated with a first-order Euler scheme with no reference to $s$ , so I'm not sure that this approach conserves the Nose energy (equation 1) -- it would be a good test to see how the scheme is performing, and if this could be the source of some stability issues. This would require keeping track of $s$, throughout the simulation.

Also this scheme gives phase trajectories that are rescaled in time, dynamic properties will be skewed, so this should be explicitely stated in the documentation and/or $s$ should again be tracked to allow for correction.

As a side note, I don't think NoseHoover needs to be a mutable struct anymore?

All the best,

Noé

[ejmeitz]

I’ll probably just add s then, and double check with y’all again. Mind if I ask you questions if I’m stuck? You definitely understand this stuff on a level I dont.

And yeah it does not need to be mutable, forgot to remove that!

[noeblassel]

Sure!

[jgreener64]

Thanks for the input Noé and for the further details Ethan.

If you can find an analogous implementation in another software that might make it easier. I've been basing a lot of the algorithms in Molly on OpenMM, for instance.

[ejmeitz]

Ok so I have not changed anything in the code, but I had the chance to read that paper closer. The equations that contain the extra degree of freedom (s) are in Nose & Hoover's original formulation, but this paper states that that DOF is unnecessary:

"We should also point out that both Nose and Hoover transformed P_s to p' = p/s. However as we shall see this additional transformation is an unnecessary complexity."

so I do not think that equation is necessary to attain a symplectic integrator. In fact I think the point of this paper is proving that this "virtual variable" is unnecessary.

I can try and prove my integrator is symplectic if needed.

Per the paper the necessary equations to integrate are:

[jgreener64]

Okay, do you think any changes are required to address this comment @ejmeitz?

except for the final velocity update which I imagine should instead be sys.velocities =v_half + (sim.dt/2) * (remove_molar.(accels_t) .- (zeta.*sys.velocities))) -- although this is consistent to the first order in sim.dt with the one proposed.

[ejmeitz]

I'm assuming that should be accels_t_dt not just accels_t. The numerical scheme I used is the exact same (leapfrog-like) scheme used in velocity verlet. Note I have not explicitly worked through the numerical scheme in awhile so I might be missing something, but I'm fairly certain my scheme is the same order as everything in Molly currently.

Heres' a source that integrates the ODE's very similar to the current code. I will double check my scheme incorporates the zeta half step in this source though.

[ejmeitz]

Adding the zeta half step made the temperature fluctuations bigger initially (see graph above which only goes up to ~340K). Not the most general test, but I dont think the half step is necessary since zeta is already being defined using the temperature at t + dt/2.

[jgreener64]

It seems more consistent to me to add the zeta half step, since zeta is used to update the velocities where the half step is used.

Then I think this is ready to merge.

[ejmeitz]

Added half step! Can add docs about it if needed.

Also my crystal package is almost done just need to get it added to a registry. When that happens I still want to change the tests for this simulator.

[jgreener64]

The zeta half steps should be divided by 2 I think as in https://www2.ph.ed.ac.uk/~dmarendu/MVP/MVP03.pdf.

Also my crystal package is almost done just need to get it added to a registry. When that happens I still want to change the tests for this simulator.

Sounds good.

[ejmeitz]

Yep, that helped the fluctuations in the protein example.