Constraints with Langevin Middle Integrator

[Ethan Meitz]

Hello,

I was reading through the integrators section and was curious to learn more about how the constraints are implemented in the LangevinMiddleIntegrator (as I am trying to code it myself). From my understanding algorithms like SHAKE and RATTLE add "constraint forces" to the existing inter atomic forces so that when you perform the velocity and then the position update the constrained atoms end up in the correct position. In the LangevinMiddleIntegrator there are two position updates and hence two applications of the constraint algorithm, but the system forces are only applied in the first update step. In the second step the system forces are not directly used rather the random forces necessary to get the Langevin dynamics are used. However, it is unclear to me how then you apply constraint forces for the second step as the system forces are not used in the second position update. Do your constraint algorithms directly act on the coordinates instead of modifying the forces acting on the atoms?

Thanks!

[Peter Eastman]

Almost all integrators apply constraints at the level of positions, and possibly also velocities. You don't integrate the constraint forces. Instead you integrate as if there were no constraints, then modify the positions to satisfy the constraints. Based on how much they change, you can back-calculate what the constraint forces must have been to produce that change in positions.

The API documentation for CustomIntegrator gives code for implementing the algorithm used by LangevinMiddleIntegrator with a CustomIntegrator.

Code: Select all

CustomIntegrator integrator(dt);
integrator.addGlobalVariable("a", exp(-friction*dt));
integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)));
integrator.addGlobalVariable("kT", kB*temperature);
integrator.addPerDofVariable("x1", 0);
integrator.addUpdateContextState();
integrator.addComputePerDof("v", "v + dt*f/m");
integrator.addConstrainVelocities();
integrator.addComputePerDof("x", "x + 0.5*dt*v");
integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
integrator.addComputePerDof("x", "x + 0.5*dt*v");
integrator.addComputePerDof("x1", "x");
integrator.addConstrainPositions();
integrator.addComputePerDof("v", "v + (x-x1)/dt");

[Ethan Meitz]

Thanks, that makes sense! So if I were to dump forces would it not include the constraint forces?

[Peter Eastman]

Correct, it would not.

[Ethan Meitz]

Another quick question on constraints.

Which algorithms are used internally by OpenMM when constraining a simulation? Reading the docs it sounds like you use SETTLE for water, but otherwise does it depend on the integrator chosen? Like for LangevinMiddleIntegrator is only SHAKE used whereas for VelocityVerlet is the constraint algorithm automatically upgraded to something like RATTLE? Thanks!

[Peter Eastman]

It depends on the platform, not on the integrator. All platforms use SETTLE for water. The Reference and CPU platforms use CCMA for everything else. The CUDA and OpenCL platforms use SHAKE for small clusters composed of hydrogens bound to a single heavy atom, and CCMA for anything else.

[Ethan Meitz]

So you never satisfy the velocity constraints in a simulation (unless its water)? Does that not create artifacts?

[Peter Eastman]

Those algorithms are just different methods of solving the same system of equations. They all produce identical results (up to the tolerance you specify). The velocities are updated correctly in all cases.

[Ethan Meitz]

Well we're far beyond my original question, so now I'm just trying to understand. The systems of equations for RATTLE and SHAKE are different though, and looking at the CCMA paper it looks like you only solve the equations generated by the position constraints. In integration schemes where you also integrate the velocity how would CCMA satisfy the velocity constraint ( dot( v_ij, r_ij) = 0 )? (thanks for humoring all my questions btw)

[Peter Eastman]

SHAKE is a constraint algorithm. RATTLE is an integration algorithm. They're different because they do different things!

Here's an implementation of RATTLE implemented with a CustomIntegrator.

Code: Select all

CustomIntegrator integrator(dt);
integrator.addPerDofVariable("x1", 0);
integrator.addUpdateContextState();
integrator.addComputePerDof("v", "v+0.5*dt*f/m");
integrator.addComputePerDof("x", "x+dt*v");
integrator.addComputePerDof("x1", "x");
integrator.addConstrainPositions();
integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
integrator.addConstrainVelocities();

Notice the line where it calls addConstrainPositions(). The job of that line is to adjust all the positions so they satisfy the constraints. That involves solving a system of equations, which can be done with various methods: SETTLE, SHAKE, etc. All of them produce the same result. There's a unique set of positions that satisfy the equations, and any constraint algorithm will find it (unless it fails because you're outside the algorithm's radius of convergence, but that's a different matter).

Conceptually the changes to the positions were caused by constraint forces. You don't need to know what those forces were in order to adjust the positions. Often it's easier to just solve for the positions directly. But the constraint forces will also affect the velocities, and that needs to be correctly accounted for. Notice that immediately before applying position constraints we copy the current positions x to a temporary variable x1. Then immediately after applying them we add a correction (x-x1)/dt to the velocities. That's a key piece of the RATTLE algorithm. It ensures the positions and velocities are both updated based on the same constraint forces.