Feature/modern integrators

[MartinGirard]

Description

Provides modern integrators to hoomd, mainly stochastic velocity rescaling & langevin piston methods. To simplify code, global thermostats / barostat are abstracted out so that all methods share GPU kernels.

Motivation and context

Currently available integrators in hoomd are "primitive". Implementation of new integrators is complicated by the need to write kernels for every integrator, even though typically only the rescaling factors change.

Resolves #1419

How has this been tested?

Will add unit tests (&fix NVT/NPT test) & integrators to lj_equation_of_state tests.

Change log

Will add the following methods

hoomd.md.methods.NVTStochastic
hoomd.md.methods.NPTStochastic
hoomd.md.methods.NPTStochasticPiston

Checklist:

• I have reviewed the Contributor Guidelines.
• I agree with the terms of the HOOMD-blue Contributor Agreement.
• My name is on the list of contributors (sphinx-doc/credits.rst) in the pull request source branch.

[MartinGirard]

@joaander I've started by abstracting parts out of the Nosé-Hoover NVT & Berendsen thermostats out into NVTBase. This scheme means there's a single kernel for all global NVT methods and its much easier to implement new thermostats/barostats. At the same time, it generates virtual inheritance patterns. Maybe there's a better way of doing this.

[joaander]

Thanks! This looks like a good start and nicely refactors the common elements of the NVT integration methods to better reuse code.

I think the virtual inheritance is a valid solution. It nicely keeps the CPU and GPU code compartmentalized. One alternative I see is to have only one NVTBase class that uses m_exec_conf->isCUDAEnabled() conditions to choose the code path. IntegratorTwoStep::update uses this method to handle computeNetForce and computeNetForceGPU. In the case of the new NVTBase, there are several methods and autotuners that differ between CPU and GPU solutions - so there will be a number of conditions that add code complexity. On balance, I think the virtual inheritance method might be the cleaner approach. Any thoughts @b-butler ?

[b-butler]

I think that trying the cleaner virtual inheritance approach at least at first makes more sense, and only if performance needs to be improved should we consider switching.

[joaander]

I am not concerned about performance. These methods are called once per timestep.

[MartinGirard]

@joaander I've added the Bussi thermostat for NVT (trivial to implement here), and started doing the same work for the NPT integrators. I initially set out to create an (MTTK-Langevin piston) +(MTTK-Langevin-barostat + Bussi thermostat) classes in addition to the base MTTK. There's already a friction factor in the base MTTK class (m_gamma). I've modified the base MTTK barostat into a Langevin piston by using the same friction term and adding a random force on the piston (setting gamma = 0 recovers MTTK). Is Langevin piston by default what we want here?

[joaander]

@joaander I've added the Bussi thermostat for NVT (trivial to implement here), and started doing the same work for the NPT integrators. I initially set out to create an (MTTK-Langevin piston) +(MTTK-Langevin-barostat + Bussi thermostat) classes in addition to the base MTTK. There's already a friction factor in the base MTTK class (m_gamma). I've modified the base MTTK barostat into a Langevin piston by using the same friction term and adding a random force on the piston (setting gamma = 0 recovers MTTK). Is Langevin piston by default what we want here?

For compatibility with existing user scripts in a v3.x release, we need to maintain the existing behavior for both gamma = 0 and gamma != 0 in the NPT and NPH classes.

We can revisit this in a v4 release in the not too distant future and decide how we want to expose the options, perhaps through the gamma parameter as you mention here or perhaps by making the NPT, NVT, and NPH base classes define only the parameter interface with specific implementations in NPTMTK, etc...

[MartinGirard]

@joaander So here's the full tentative code. I pondered for a while on which class hierarchy was appropriate for NPT. Only the thermostat changes between the different barostat, so the base class now implements pure MTTK nph integration. Implementation should provide backwards compatibility with the rest of hoomd 3.x.

For v4, I would propose deprecating the current m_gamma implementation, so that gamma != 0 automatically yields a Langevin piston. The base class could then implement NPH, with derived classes implementing thermostatting.

Code still needs more documentation / proper formatting in some places, and I have to fix the unit test (they are currently broken in hoomd's repo). I think it would be a good time for comments on the implementation before I fix docs / formatting, also class names. A particular point of interest: the Bussi thermostat divides by the current temperature to compute rescaling. There's no verification that the user did any thermalization, and if velocities are zero, then this will result in infinite rescaling factors. Do we want to provide any check there?

I have a side question. I had to add zero-initialization for the barostat (m_barostat struct) in the NPT class. When looking at the current hoomd repo, I can't find any place where its properly initialized, and this confuses me a bit.

[joaander]

@joaander So here's the full tentative code. I pondered for a while on which class hierarchy was appropriate for NPT. Only the thermostat changes between the different barostat, so the base class now implements pure MTTK nph integration. Implementation should provide backwards compatibility with the rest of hoomd 3.x.

For v4, I would propose deprecating the current m_gamma implementation, so that gamma != 0 automatically yields a Langevin piston. The base class could then implement NPH, with derived classes implementing thermostatting.

I will review the code soon and offer comments.

After we merge this pull request, we can make the v4 breaking changes via another pull request off the trunk-major branch. I expect a v4 release in ~6 months: https://github.com/glotzerlab/hoomd-blue/issues?q=is%3Aissue+is%3Aopen+label%3Abreaking

Code still needs more documentation / proper formatting in some places, and I have to fix the unit test (they are currently broken in hoomd's repo). I think it would be a good time for comments on the implementation before I fix docs / formatting, also class names.

Are you referring to the C++ unit tests? Do you think they are worth keeping? In v3, we test basic functionalities in Python via pytest. I'm currently working on a separate repository with longer running validation tests that ensure simulation results are statistically correct. I will build this branch and test it there.

A particular point of interest: the Bussi thermostat divides by the current temperature to compute rescaling. There's no verification that the user did any thermalization, and if velocities are zero, then this will result in infinite rescaling factors. Do we want to provide any check there?

Yes, we should provide the user with a meaningful error message via an exception.

I have a side question. I had to add zero-initialization for the barostat (m_barostat struct) in the NPT class. When looking at the current hoomd repo, I can't find any place where its properly initialized, and this confuses me a bit.

Adding zero-initialization in C++ is good practice, keep it. The current HOOMD release does initialize m_barostat from Python, but the code path is not obvious. barostat_dof is in the NPT classes param_dict. When _HOOMDBaseObject initializes the C++ class, it calls self._apply_param_dict() which in turn will call setBarostatDOF with the values provided by the user, or the default of (0,0,0,0,0,0) of the user didn't provide a value.

[MartinGirard]

Yes, we should provide the user with a meaningful error message via an exception.

I've added an exception there, as well as in the Berendsen thermostat.

Are you referring to the C++ unit tests? Do you think they are worth keeping? In v3, we test basic functionalities in Python via pytest. I'm currently working on a separate repository with longer running validation tests that ensure simulation results are statistically correct. I will build this branch and test it there.

I wasn't aware that you were deprecating the unit tests. Its likely better to guarantee that the algorithm produce correct ensembles.

I've also made that one header consistent with virtual / override keywords. Should the code base favor override?

[joaander]

Are you referring to the C++ unit tests? Do you think they are worth keeping? In v3, we test basic functionalities in Python via pytest. I'm currently working on a separate repository with longer running validation tests that ensure simulation results are statistically correct. I will build this branch and test it there.

I wasn't aware that you were deprecating the unit tests. Its likely better to guarantee that the algorithm produce correct ensembles.

Most of the C++ unit tests were originally written before HOOMD had a Python interface. I'm slowly replacing most of them with tests written in Python. There are still some things that need to be tested in C++, but most are best tested in Python. Specifically, anything that requires initialing a full system state and running a small number of timesteps are best done in Python. The full policy is defined here: https://hoomd-blue.readthedocs.io/en/v3.6.0/testing.html#implementing-tests

The ensemble testing code isn't available publicly yet, but it will be soon.

[joaander]

pre-commit.ci run

[joaander]

pre-commit.ci autofix

[joaander]

Thanks for all your work on this!

I added documentation comments to Thermostat, TwoStepConstantVolume, and TwoStepConstantPressure.

After fixing the rigid body bug in #1469, the rigid body integrators in #1424 pass the physical validation tests.

[b-butler]

Thanks @MartinGirard for all the work. I left one correction (I think) to a comment but this is ready to merge in my opinion.