Custom integrator for diamondoid machinery

[philipturner]

I am trying to evaluate bonded and nonbonded forces at indivisible frequencies. With MM4, the stretch, bend and bend-bend (hereafter "bonded") are stable at 2.5 fs. Meanwhile, the nonbonded, nonbonded 1-4, torsion, torsion-stretch, and bend-torsion-bend are stable at 4.0 fs. I wish to use a staggered 2:3 sample rate, running the bonded forces at 2.67 fs instead of the current 2.0 fs. However, I am at a loss for how to debug a symplectic integrator for this.

I tried this and the energy instantly became NAN. I am not using any constraints or a thermostat (yet), because the existing ones violate conservation of angular momentum in localized areas.

integrator.addGlobalVariable(name: "initialized", initialValue: 0)
integrator.addPerDofVariable(name: "nonbonded", initialValue: 0)
integrator.addPerDofVariable(name: "bonded", initialValue: 0)

let loopIterations = Int(exactly: Self.timeStepInFs / 8.00)!
for _ in 0..<loopIterations {
  // Create a time-reversible sequence of integration steps.
  var sequence: [() -> Void] = []

  integrator.beginIfBlock(condition: "initialized < 1")
  integrator.addComputeGlobal(variable: "initialized", expression: "1")
  integrator.addComputePerDof(variable: "nonbonded", expression: "f1")
  integrator.addComputePerDof(variable: "bonded", expression: "f2")
  integrator.endBlock()

  // Bonded force evaluation (1 / 3).
  sequence.append {
    integrator.addComputePerDof(variable: "v", expression: """
      v + (1 / 6) * (dt / \(loopIterations)) * bonded / m
      """)
  }
  sequence.append {
    integrator.addComputePerDof(variable: "x", expression: """
      x + (1 / 6) * (dt / \(loopIterations)) * v
      """)
  }
  sequence.append {
    integrator.addComputePerDof(variable: "bonded", expression: "f2")
  }

  // Nonbonded force evaluation (1 / 2).
  sequence.append {
    integrator.addComputePerDof(variable: "v", expression: """
    v + (1 / 4) * (dt / \(loopIterations)) * nonbonded / m
    """)
  }
  sequence.append {
    integrator.addComputePerDof(variable: "x", expression: """
      x + (1 / 12) * (dt / \(loopIterations)) * v
      """)
  }
  sequence.append {
    integrator.addComputePerDof(variable: "nonbonded", expression: "f1")
  }

  // Prepare for bonded force evaluation (2 / 3).
  sequence.append {
    integrator.addComputePerDof(variable: "v", expression: """
      v + (1 / 3) * (dt / \(loopIterations)) * bonded / m
      """)
  }
  sequence.append {
    integrator.addComputePerDof(variable: "x", expression: """
      x + (1 / 4) * (dt / \(loopIterations)) * v
      """)
  }
  for statement in sequence {
    statement()
  }

  // Center of time-reversible sequence.
  integrator.addComputePerDof(variable: "bonded", expression: "f2")
  integrator.addComputePerDof(variable: "v", expression: """
    v + (1 / 2) * (dt / \(loopIterations)) * nonbonded / m
    """)

  for statement in sequence.reversed() {
    statement()
  }
}

[peastman]

The rRESPA algorithm requires the outer step to be an integer multiple of the inner step. In principle it should be possible to come up with a different way of factorizing the propagator that allows overlapping steps like you described. The theory is described in https://doi.org/10.1063/1.463137. I'm doubtful that it would work, though. You would effectively have a compound step with an outermost step size of 8 fs. MTS integrators tend to have resonances that make them unstable unless the outermost step is less than half the period of the fastest motion in the system. In typical molecular simulations with constraints on hydrogens, that limits them to about 4 fs.

[philipturner]

I was able to fix a bottleneck in the bend-bend force, so the cost of bonded forces is less significant. I settled on the following scheme:

• bonded forces: 2.17 fs (100 / 23 / 2)
• nonbonded and torsion forces: 4.35 fs (100 / 23)

What’s your recommendation for a thermostat that would conserve angular momentum of a rigid body?

[peastman]

A thermostat models the interaction with a heat bath. If the system is in contact with a heat bath, its angular momentum shouldn't be conserved. Can you explain what you're trying to do?

[philipturner]

I want to run nanosecond-long simulations of simple nanomechanical systems, around 10K-100K atoms. The issue is, energy starts to explode after a few hundred picoseconds. I need a way to make it have reasonable kinematic behavior at these time scales, without resorting to smaller time steps.

Before the simulation starts

To start, I am calling setVelocitiesToTemperature before the simulation starts. I retrieve and correct the velocity of each rigid body individually (via OpenMM_State), to set its atoms' net linear and angular momentum to zero. Then I add a specific (angular) velocity that I actually want it to have.

During the simulation

At this stage, rotating parts are not manufacturable yet, so conserving each part's angular momentum w.r.t. other parts may not be a priority (design the system to remove rotational degrees of freedom). However, I don't want the entire system spinning while I try to visualize it. An absolute necessity is relative linear momenta. That should be easier to conserve in the thermostat (constraints maybe?).

[peastman]

CMMotionRemover computes the total linear momentum of the system at each step and sets it to zero. You could create something similar for rotational motion. At each step, compute the total angular momentum around some point and set it to zero. That won't completely prevent the system from rotating, since you can still get zero angular velocity rotations, but those should at least be much slower.

Before doing that, though, I would verify that you can get good energy conservation. If you use smaller steps, mixed precision mode, and a tight tolerance for constraints, you should have no trouble simulating hundreds of nanoseconds with minimal energy drift. If you can't, that suggests something else may be wrong.