Angular momentum subtraction in Langevin dynamics

[rodrigo]

Hi,

I wrote a updater plugin to zero out the angular momentum and angular acceleration of a collection of confined point particles that are being integrated using Langevin dynamics and I am wondering if I may be missing something subtle.  The plugin seems to work reasonably well, but, depending on when the update is executed, zeroing out the angular acceleration may or may not make any difference.  I guess my concrete question is then: when are updaters executed?  As I understand it, accelerations are computed from scratch during second integration step, so I would lose any corrections that I compute to them if the updaters were run between the first and second steps: which may fine for updaters that act only on particle velocities.  However, if the updater executes after accelerations are computed in step two, then I should be able to maintain a strictly irrotational system for all time.

Also, presumably the effective temperature will now be reduced to T*(1-2/ndof), but that's easy to fix.  Is there anything potentially subtle that I should worry about?

BTW, it seems to me like ZeroMomentumUpdater should be augmented to zero out forces as well.  Newtonian integrators obviously satisfy this, but the bd forces in the langevin integrator will push the center of mass around with a displacement proportional to deltaT^2.

Best,

Rodrigo

[Joshua Anderson]

Evaluation order is:

* Analyzer1

* Analyzer2
* ...

* AnalyzerN

* Updater1

* Updater2
* ...

* Updater N

* Integrator

    * Step one

    * Compute force1, force2, ... forceN

    * Compute net force

    * Compute constraint_force1, constraint_force2, ... constraint_forceN

    * Modify net force due to constraint forces

    * Step two (may involve computing accelleration)

With this order, analyzers and updaters always work on x(t) and v(t) and do not need to be concerned about half steps for any quantities.

The acceleration particle data field is only cached from the end of step two to the beginning of step one. If you modify it in an updater, you only modify the acceleration applied in step 1, not step 2. Essentially, the Integrator owns acceleration and should be the only class to touch it. To modify accelerations in both steps, you need to implement a ForceConstraint which runs with knowledge of the forces on particles due to non-constraint forces. Constraint forces advertise the number of degrees of freedom that they remove for proper temperature accounting.

For your specific request on the Langevin integrator, a ForceConstraint is not sufficient. Langevin implements its random force within the integrator itself, it is not a ForceCompute. I have long considered adding an option that zeroes the "net random force" on particles due to the langevin thermostat each step, but have not taken the time to do it. The correct place to implement that is within the integrator itself. Integration methods provide getNDOF which advertises how many DOFs they add to the system for proper temperature accounting. Pull requests are welcome.

As for ZeroMomentumUpdater, it predates the langevin thermostat and is vestigal anyways. I implemented it to remove momentum drift from single precision NVT/NVE integration on GPUs with poorly rounding FMADs. Now that GPUs all have FMA and double precision, NVT/NVE drift is practically undetectable. (In single precision, it is possible to detect, but only in extremely long runs).

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To post to this group, send email to hoomd...@googlegroups.com.
Visit this group at https://groups.google.com/group/hoomd-users.
For more options, visit https://groups.google.com/d/optout.

[Rodrigo Guerra]

Thanks, that is very helpful.  I should probably rewrite my updater as a constraint; but, since the Langevin integrator does compute accelerations, this works.

Best,
Rodrigo

[Åsmund Ervik]

Just a comment: ZeroMomentumUpdater (as a one-shot thing) is very useful when one applies an external periodic potential e.g. to speed up phase separation. Since the phases (different particle types) will typically have a center-of-mass not aligned with the center of the external potential well, the system can gain a net momentum which makes post-processing much more difficult. Being able to remove the net momentum is very handy.

On Thursday, February 4, 2016 at 2:11:31 PM UTC+1, Joshua A. Anderson wrote:

 ...

[Rodrigo Guerra]

Quick follow-up to this:  Since the random forces of the BD and Langevin integrators are computed in step two, I imagine that they will violate the constraints.  At least for ConstraintSphere, the positions are self-correcting, but the effect on the velocities is O(sqrt(deltaT)), which probably affects the reported KE.  In this case it may make sense to apply the constraint as an updater that fiddles with velocities and accelerations.

Does that seem right?

Thanks,

Rodrigo


On 02/04/2016 08:11 AM, Joshua Anderson wrote:

[Joshua Anderson]

No, it does not. Even if it were possible to insert updates between step1 and step2, the pressure computation would be off as the constraint force would not be properly accounted for.

To properly account for constraint forces such as ConstraintSphere with Langevin and BD dynamics, the constraint method needs to be taught how to predict what the integrator is going to do on the particle. The constraint forces work like this: Given the net force before the constraint F_n, and the knowledge of how the positions are updated:  x(t+1) = f(x(t), F_n), solve for a force F_n' that so that x'(t+1) = f(x(t), F_n') satisfies the constraint. To do this, it needs to pretend to update a step, project x down onto the constraint, then solve for the additional constraint force F_c = F_n' - F_n that will cause the constraint to be upheld.

Currently, ConstraintSphere assumes Velocity Verlet dynamics (integrate.nve). Given a python-level link to the integrator, it could retrieve the gamma values, seed, etc... needed to determine the appropriate constraint force for the other integrator types. It is just a matter of solving the appropriate equations. Random numbers are not a problem because of the hash based RNG.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

[Rodrigo Guerra]

I see.  Thanks for the clarification.  I chose to mitigate this problem to first order with a clone of TwoStepLangevin that zeros out dot(pos, bd_force), which reduces the error in the pressure and variance of KE  dramatically; however, as you point out, this correction is not exact and a proliferation of ancillary updaters/integrator clones is worse than fixing the problem as you suggest.

Best,
Rodrigo

[Joshua Anderson]

On revisiting this, I realize now that the plan I put forward is not a workable solution. Langevin/brownian dynamics with constraints needs to be implemented with the knowledge of the constraint in the integrator. Case in point, with the sphere constraint, the particle is on a 2D surface embedded in 3D. The random force should be generated in 2D perpendicular to the surface, but the code currently computes a 3D random force which is effectively projected down onto the surface by the constraint force. This modifies the distribution of random forces, resulting in different temperature than the set point.

I made an issue to track this: https://bitbucket.org/glotzer/hoomd-blue/issues/132 . I hope to include it in release v2.0.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

[Rodrigo Guerra]

Hi Joshua,
I had not appreciated this earlier, but you're absolutely correct for a force distribution uniformly distributed in a cube; however, there's no problem with projection if the distribution is isotropic.  I imagine that the easiest solution is to keep the RNG, but to generate numbers until the force satisfies rx*rx+ry*ry+rz*rz<=1, which should still be deterministic and happen a little over 50% of the time.
-R

[Rodrigo Guerra]

Hi Joshua,
Come to think of it, the rejection method is probably pretty lousy for GPUs if all the threads have to wait for the one laggard that happens to take 20 tries to find an adequate vector.

If you want to be slick, you can generate 6 random numbers on [-1,1] and multiply each of them by one of six icosahedral vertex directions (12 modulo inversions).  This generates an isotropic covariance matrix with a displacement distribution that looks a lot like a gaussian.

Best,
Rodrigo

On 04/14/2016 07:04 AM, Joshua Anderson wrote:

[Rodrigo Guerra]

Hi Joshua,
I should amend my previous comments.  The variance of the Langevin force projected along any arbitrary axis is (surprisingly?) isotropic for either the uniform cube or the icosahedral sum.  However, the kurtosis of the projected force is strongly anisotropic and quite negative for the cube, whereas it is rigorously isotropic and about half as large for the icosahedral sum.

That's to say, no matter how the constraint algorithm projects the stochastic force, the average temperature will be fine: which is consistent with all the thermodynamic results I have on the sphere.  Moreover, if dt is small, I would expect that the anisotropy in the kurtosis won't be a problem.  Then again, implementing the icosahedral sum is so easy that it's probably worth the extra computation.

Best,
Rodrigo


On 04/14/2016 07:04 AM, Joshua Anderson wrote:

[Isaac Bruss]

I believe the issue isn't with generating a distribution uniformly within a sphere rather than a cube, but is instead the problem of projecting the distribution, even if it is correctly uniform within a sphere, onto a surface. For instance, a random force that happens to be normal to the surface, when projected onto that surface will yield a force of zero. So while the distribution is ideal in 3D, it has "shrunk" after being projected, and becomes narrower than desired.

peace,

- isaac

[Rodrigo Guerra]

Sure, but a force that points straight up has measure zero. What's important is that the variance of the projected force satisfies the 2D FDT: which is the case even for forces sampled uniformly from a cube and projected onto an arbitrary axis or plane. The 1D projection of a vector randomly sampled from inside a sphere is obviously isotropic, it's PDF is pretty close to cos(pi F/2), has finite moments, and all that.

The advantage of the icosahedral sum is that the variance and fourth moments are both isotropic, and the anisotropy of the sixth moment is tiny (<1%). Its PDF also looks pretty Gaussian, albeit with finite support. So you can be fairly certain that the orientation of the constraint surface doesn't matter. Then again, I've never seen a problem with the thermodynamics (average and variance of KE or specific heat) of particles on a sphere with the current random force. But, I could imagine that using very large time steps might not give the central limit theorem enough iterations to wipe the anisotropy of the higher moments in some time scale relevant to the dynamics: which is really user error.

R

[Joshua Anderson]

I am still concerned because the force is not actually projected to the 2D surface, but a 3D move is made to the particle based on the equations of motion, and then the particle is projected back to the surface by some means. This is only effectively projecting the random force. I think it would be more robust if the random force was generated in the proper dimensionality first and doing so is not a significant amount of development effort.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

[Isaac Bruss]

I see your point, Josh. Plus, I think my previous interpretation was wrong. A particle confined to a surface doesn't actually live in 2D, but rather is still in 3D with a strong adhesion to the surface. Therefore the distribution shouldn't be pure 2D, but rather 3D and then projected onto the surface. And as you are noticing, this projection step isn't the most robust in its current implementation.

- i

--
You received this message because you are subscribed to a topic in the Google Groups "hoomd-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/hoomd-users/yrSAr60pJsQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to hoomd-users...@googlegroups.com.

[Rodrigo Guerra]

Fair enough. My kludgey implementation still uses constrain.sphere, but modifies step two of the integrator so that it generates a 3d random force, projects out the normal component, and adds a tiny inward force to adjust the centripetal force. This is very simple and works very well. I'd be happy to send you the relevant plugins if that's helpful at all.

-R

[Rodrigo Guerra]

I don't quite see why the langevin force is more problematic than the net normal force due to deterministic pair interactions.  As I understand constraint force class, there's nothing that guarantees that the force required to project the unconstrained trajectory onto the closest surface point will be orthogonal to the displacement: which is necessary for conservation of energy.  For example, the constraint force computed by the nearest point projection method on the sphere will be proportional to the surface normal at the end point; if you have a system composed of strongly repulsive particles, as I do, the constraint force will be large and negative and will do work against the motion of the particle.

You can solve this to leading order by setting the normal force to zero before handing the problem over to constraint.evalU(), and you can do better if you know the local curvature tensor.  So, langevin force or not, the net effect of this would be to project 3D forces onto the 2D tangent plane.

-R