In short, I am trying to implement my own constraint to NPT dynamics in ASE (fixing center of mass, as well as no-overall-rotation. The latter is not implemented in the existing constrain classes, I think). My naive definition of my own constraint class didn’t work, so I’m trying to understand the how constraint is applied in NPT source code. However, I couldn’t find the places where the constraint is implement in the NPT source code. Please tell me which line the constrains are implemented.

In details, I can identify how other molecular dynamic handle constrains.

For instance, I found VelocityVerlet enforce constrain in the line “atoms.set_positions(r + self.dt * p / masses)”. Langevin has “atoms.set_positions(x + self.dt * self.v + self.rnd_pos)” where the set_positions from the atoms class will enforce the constrains.

The closest thing I found in NPT that apply constraint is

"

def _setbox_and_positions(self, h, q):

“”“Set the computational box and the positions.”“”

self.atoms.set_cell(h)

r = np.dot(q + 0.5, h)

self.atoms.set_positions(r)

"

but the r = np.dot(q + 0.5, h) doesn’t make sense to me. It doesn’t look like the expect q=q_old+ t * v.

I’m not very familiar with Nose-Hoover thermostats. Is the above code where the constrain is implemented?

Thanks,

Bohan