[WIP] Added bitwise-reversible velocity Verlet integrator

[jchodera]

I've added a draft of a bitwise-reversible velocity Verlet integrator along with a battery of tests to make sure that the positions end up back where they started after 10 steps of integration, velocity reversal, and 10 steps of integration.

Using the integrator looks like this:

# Create a bitwise-reversible velocity Verlet integrator.
integrator = integrators.BitwiseReversibleVelocityVerletIntegrator(timestep)
# Create a context.
 context = openmm.Context(testsystem.system, integrator, platform)
# Set positions and velocities.
context.setPositions(positions) 
context.setVelocities(velocities) 
# Truncate accuracy of positions and velocities.
# This is required for bitwise-reversibility!                                                                                                                                            
integrator.truncatePrecision(context)
# Integrate dynamics in a bitwise-reversible manner.
integrator.step(nsteps)                                                                                                                                      

There is control over two kinds of precision:

# Create a bitwise-reversible velocity Verlet integrator.
integrator = integrators.BitwiseReversibleVelocityVerletIntegrator(timestep, precision_nbits=20, timestep_precision_nbits=5)

The quantity precision_nbits is the number of bits of precision retained in the quantity 0.5*dt*f/m accumulated in the velocity of velocity Verlet. This is done crudely via the truncation

floor(0.5*dt*f/m*scale + 0.5)/scale

where scale = 2**precision_nbits, which means that quantities smaller than 2**-precision_nbits in natural OpenMM velocity units (9.5e-7 nm/ps in this example) are lost. Note that this transformation does not really retain precision_nbits bits in the mantissa, and is more appropriate to fixed-point or integer arithmetic.

The quantity timestep_precision_nbits is used to round the specified timestep to retain only this number of bits in the floating-point mantissa. If you want a timestep of 1 fs, for example, the best representation with 5 bits of mantissa precision (~0.000977 ps) will be retained. This is done so that the position update step of velocity Verlet

x+dt*v

will not accumulate too many bits of precision during integration.

Note that this scheme is not ideal, as I fear it is possible to accumulate timestep_precision_nbits bits of precision in the positions and velocities each timestep. Even with double precision accumulation, this could quickly overwhelm the number of available bits of precision. I'm working on finding a way to rework the velocity-explicit velocity Verlet scheme to avoid this issue.

There are some limitations on what kinds of systems and platforms are supported:

• Only platforms with deterministic forces (Reference, OpenCL) are supported
• Only precision models in which positions and velocities are accumulated in double precision (mixed and double) are supported
• Constraints are not supported; handling constraints in a bitwise-reversible manner seems to be an open problem for the field

[jchodera]

There is more information on how to implement this correctly in a patent filing from DESRES:
https://www.google.com/patents/US7809535

I'll digest this and take another stab at this soon.

[jchodera]

I tried to implement the scheme from the DESRES patent, but the rounding applied to the position kick in this line doesn't seem to give me bitwise-reversibility. Only my scheme with the precision-degraded timestep and no rounding on the position update seems to work, though I am still worried that this scheme will continue to accumulate significant bits of precision each timestep and quickly overwhelm even double precision position and velocity accumulators.

I suspect this might be due to the use of floor(scale*s + 0.5)/scale to degrade the precision of quantity s, since this may not be exactly the same as the rounding operation required (where rounding following scaling should be done to the nearest integer with ties broken by rounding to the even integer), but I am not certain of this.

[kyleabeauchamp]

I am not a lawyer, but shouldn't we be avoiding this patent with a 10 foot pole? I thought even reading the patent puts us at risk of a civil lawsuit where they claim we are using patented technology without a license.

Our lawyer friend @rmcgibbo might have some insight on the matter. The extent of my knowledge is from quora: http://www.quora.com/Is-it-illegal-to-implement-a-patented-algorithm

[kyleabeauchamp]

I realize that the "clean room" stuff is more for copyright infringement, but still AFAIK we can't implement and distribute their patented algorithm.

[jchodera]

• I doubt this scheme was actually patentable.
• We are not making any money off of its implementation
• I really honestly don't care

[jchodera]

Amusingly, the patent contains so many typos that it nowhere describes a valid scheme. I wonder if that was intentional. Technically speaking, we cannot actually implement the scheme described in the patent.

[kyleabeauchamp]

Sounds good. Sometimes it's better to be safe than sorry with these things, if only as a way to avoid having to deal with legal headaches / timesinks.

[jchodera]

Fair enough. To be clear:

• What is described in the patent is not what is implemented here
• I experimented with implementing a concept hinted at by the patent in which the position update was also subjected to fixed-precision rounding, but that does not appear to produce bitwise-reversible trajectories
• The bitwise-reversible velocity Verlet scheme described in equations in the patent is not a valid algorithm, hence we could never implement it
• The correct scheme is not fully described by other means, therefore I don't think the patent could apply to any specific scheme we implement unless someone tries to argue that this covers all "Large floating point code, scalar or parallel, that ensures exact reversibility for non-trivial dynamics simulation", which I doubt is possible without describing a correct algorithm even if some of the ideas are right. But I am not a lawyer.
• The whole idea of trying to patent this is kind of bullshit.

[jchodera]

Maybe we should patent the correct algorithm so they need to license it from us to run Desmond...

[mrshirts]

• The whole idea of trying to patent this is kind of bullshit.

Yeah, gotta agree with you on this. Math (in theory) can't be patented.
An integrator is by definition math.

The issue is the amount of time and money required in order to fight the
B.S.

I would argue that in most cases, a bitwise reversible integrator isn't
worth it.

[jchodera]

I would argue that in most cases, a bitwise reversible integrator isn't worth it.

My suspicion is that it is necessary to make Metropolized integration work correctly for these single-precision GPU-based force calculations. Just working toward testing that hypothesis right now.