Velocity dependent forces in hoomd

[lassejan]

Hi, I have a question regarding a custom force, which is dependent on particle velocities. If I were to use the hoomd.md.methods.ConstantVolume()-integrator, which as far as I understood it uses a Velocity-Verlet integration scheme, to evolve translational degrees of freedom in time, at what time step would velocities be evaluated and used for force updates? Thank you in advance for your answers.

[joaander]

In Velocity-Verlet style integrators, including the one implemented in HOOMD-blue, forces are computed after the first half step and before the second.

[Beuteltee]

Hey i also have a bit of trouble with velocity dependent forces using the ConstantVolume integrator.
I have a driving force defined as $\textbf{F}=-\textbf{v}\frac{k}{s}(|\textbf{v}|-s)$ where $k$ is the strength and $s$ is the target velocity for the particle. I also have implemented a force cap which holds the absolute of the force vector at $|\textbf{F}=800$ to avoid crashes.
The simulation however still crashes for small values of $s= 1e-5$ with big values of $k=2000.0$. The initialisation is on a grid with random uniformly distributed velocities all normalized to $|v|=1$ so i understand for this configuration my system has a lot of energy.
What i see is that in the first step from my init. conditions to the first integrated frame the system gains a lot of velocity, so high that the simulation crashes instantly. I checked this with a smaller integration step where the simulation was stable and also saw there that the velocities are getting smaller until the system freezes completly. This is what i was expecting because i can observe this dynamic already for smaller ratios of $k/s$ with small but bigger $s\approx 0.1$. So i am curious why this is happening, since it seems to me that the integrator works fine, but after or at the initialisation i get into troubles.
Thank you in advance for your answer.

[joaander]

HOOMD-blue makes no guarantees that the equations of motion that it attempts to integrate are stable. Apply too large a force or (equivalently) use too large a value for dt and you should expect the numerically integrated solution to diverge exponentially away from the real solution. Typically, this results in a particle out of the box error in HOOMD-blue.

Making this problem worse is that your force is velocity-dependent. The derivation of Velocity-Verlet assumes that force depends only on position. By computing a velocity dependent force, you violate this assumption and introduce yet more error in the solution. Many MD codes, including HOOMD-blue, allow velocity dependent forces because those are typically small and the errors absorbed into the thermostat.

If you have need to integrate such a stiff equation of motion, I suggest you implement a custom code and choose an appropriate integrator, possibly one that uses a variable time step. Most MD codes implement only symplectic integrators due to their many nice properties, including excellent energy conservation, and typical MD simulation use-cases are well-behaved. There are other solutions (such as energy minimization or MC simulation) to resolve large position-dependent forces in the initial condition.