Fix the center of mass in MD simulation

[thomasloux]

It would be great to have a fix center of mass in MD, by removing total momentum during the simulation. It is the case at initialization but it may shift for Langevin during the simulation. I may add it later, but I wanted to make sure you haven't implemented it yet.

[curtischong]

I do not believe that anyone has implemented this!

[thomasloux]

Ok I'm trying to find on other implementations how to do that carefully. Most of the implementations fix the center of mass by fixing the center of mass velocity. This is actually even more important, with issues like flying ice cube (kinetic energy transfered from from internal system to external: https://en.wikipedia.org/wiki/Flying_ice_cube) and you would still suffer from this even by shifting then center of mass position. This also affects the computation of the temperature, based on the equipartition theorem (<E_kin> = N_{degrees of freedom)/2 * kB*T.

So looking at Gromacs (https://manual.gromacs.org/documentation/2021.1/reference-manual/algorithms/molecular-dynamics.html#kinetic-energy-and-temperature) and LAMMPS (https://docs.lammps.org/compute_temp.html): they both actually compute the temperature directly by assuming null center of mass velocity (N_DOF = n_dim N_atoms - n_dom - N_constraints). Notice that they can also account for other degrees of freedom from constraints, features not yet implemented in TorchSim. ASE by default does not assume that, but it can be enforced through Constrain (FixCOM: https://gitlab.com/ase/ase/-/blob/master/ase/constraints.py?ref_type=heads#L253). Notice however a weird (technical debt?) implementation for their Langevin thermostat that accepts a fixcm argument. They remove directly in the integrator and scale the velocities to correct the temperature computation (that then does not about the reduced degrees of freedom). For me this approach seems wrong but I don't have a careful argumentation against it. (ref PR: https://gitlab.com/ase/ase/-/issues/1285, https://gitlab.com/ase/ase/-/merge_requests/3567).

I suggest to add a Constraint class like ASE and change the computation of degrees of freedom and temperature from the quantities.py file (reintegrate degrees of freedom in SimState and add a constraint property?)

What do you think about that?

[orionarcher]

Hi @thomasloux. Thanks for looking at different implementations across the ecosystem. I agree we definitely need a constraints API, both for this use case and other simulation types.

I suggest to add a Constraint class like ASE and change the computation of degrees of freedom and temperature from the quantities.py file (reintegrate degrees of freedom in SimState and add a constraint property?)

Could these be implemented separately or do we need constraints before fixing the COM issue?

[thomasloux]

Hey, I think we need constraints to fix the COM, so that it works for any algorithm.
Still, regarding Langevin, I actually think the correct way to do that is to reproduce this paper: https://arxiv.org/pdf/1006.4914
See also Molecular Dynamics, Leimkuhler & Matthews for more details and bigger perspective (4.3 Methods for Constrained Dynamics and Constrained Langevin Dynamics). They give the correct algorithms for any constraint. I'll probably just implement it in the case of COM for Langevin. But then I need to make it would work correctly with the Constraint implementation.
All of this to say, I still need to think about it carefully.

[thomasloux]

To be honest, all of this is quite technical and difficult. Main pain points are to be sure to preserve symplecticity and conserve the order of the integrator, most of the ones we use are second order and I believe it's the case of the ones implemented in TorchSim.
Actually, in the case of COM, this correspond to a linear constraint, which have close form projector (directly project to the desired subspace: eg. remove COM momenta in the atoms momenta). So that I think they should be no issue.
In a more practical point of views, ASE constraints are applied on all integrators, although they clearly justify it only for Langevin by providing the reference and equations with and without constraints. Similarly, I have a reference for constrained BAOAB, splitting that is used in TorchSim Langevin (Leimkuhler B, Matthews C.2016 Efficient molecular dynamics using geodesic integration and solvent–solute splitting. Proc. R. Soc. A 472: 20160138.http://dx.doi.org/10.1098/rspa.2016.0138). I think it should be ok for Nosé Hoover as well. But I don't know for current NPT integrator.

So, I suggest transferring ASE implementation of linear constraints (fixAtoms, FixCOM..). But this implies:

• adding a constraints variable in SimState
• Implement a setter and getter for SimState (which shouldn't be hard for a dataclass)
• Reintegrate degree_of_freedom in SimState and modify any code computing DOF as n_dim*n_atoms in the repo (eg. in Nosé Hoover Chain, computation of temperature in general)
• Either use the same mechanism of adjust_positions, adjust_forces, adjust_momenta: use adjust_positions and momenta for MD and adjust_positions and forces for Structural Optimization. Or only use constraints on MD.

[thomasloux]

Fixed in #294