Strange result when using MonteCarloAnisotropicBarostat

[zhenyuwei99]

Hi, I am using MonteCarloAnisotropicBarostat to relax a carbon nanotube system in the NPT ensemble. In my system, the box is split by two pieces of graphene and connected with a carbon nanotube, between which is a vacuum. So I use the MonteCarloAnisotropicBarostat to relax the system's pressure only in the z-direction by creating a barostat instance like

pressure = 1 * unit.atmosphere
barostat = openmm.MonteCarloAnisotropicBarostat(
    openmm.Vec3(pressure, pressure, pressure),
    temperature,
    False,
    False,
    True,
    25,
)
system.addForce(barostat)

To fix the carbon nanotube and graphene in the correct positions, I also add a constraint force on all carbon particles like:

restrain_force = openmm.CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2")

And I run the simulation for 1ns, the result is strange as shown below: After simulation, the system gives some water cavities near the graphene. However, in some other systems, the cavity also appears in the box boundary like

I am wondering which part is incorrect in this simulation and whether the MonteCarloAnisotropicBarostat is suitable in this case

Best,

[peastman]

I'm not sure exactly what problem you're describing. The images show lots of water on either side of the tube, and a smaller amount of water passing through it. I assume that's what you expect? The irregular shape of the water on either side suggests it hasn't been wrapped to a single periodic box, so I can't tell much about the distribution of water within the box. It looks like the density of water right next to the graphene may be reduced, but since it's hydrophobic that's not surprising.

What should I be looking at?

[zhenyuwei99]

Sorry for the misleading description. Actually, I have already wrapped all the positions in the single pbc box. However aqueous box is not uniformly distributed, while some dilute volume exists in the result as shown in the picture denoted by the red curve

[peastman]

That looks like a bubble, which usually indicates the box is too big for the amount of water. The barostat ought to fix it. Can you monitor the box size? Is it getting smaller?

[zhenyuwei99]

It is getting smaller, the box size in z-direction shrink from 100A to 94A but it seems still to big

[zhenyuwei99]

When I keep increasing the simulation time to 25ns, the result seems better

The bubble size is smaller but I think still needs more than 10ns to converge as the box volume still shrinks.

Maybe relaxing pressure in only one direction has a big influence on the relaxing efficiency. I am wondering if there is any better choice for this condition. I notice there is a MonteCarloMembranceBarostat. But I am not sure whether is suitable in this case and how to determine the surface tension in my system.

[peastman]

That looks much better. This is just a guess, but having such a big completely rigid block right in the middle may be making things harder for the barostat. It still works, but that forces it to use very small steps.

Reducing the barostat interval so it takes more frequent steps might help. If you change it from 25 to 10, that might help it relax faster. It's probably only an issue for equilibration, though. Once you get into the production stage, 25 should be fine. And it looks like you've successfully equilibrated it.

[zhenyuwei99]

Good idea. We don't need the dynamic to be completely physical in the equilibrium part. When I turn the barostat frequency to 5, it converges faster.

Thanks a lot!

[zhenyuwei99]

It seems that decreasing the barostat frequency will also decrease the simulation speed. For the same system, the simulation speed decreased from 345ns/day to 243ns/day when changing the barostat frequency from 25 to 5

Maybe the barostat requires reading data from GPU and modifying the positions of particles?

[peastman]

Every barostat step requires two energy evaluations, which makes it roughly equivalent to two time steps in cost. That's why you don't want to do them too frequently. Once you're done with equilibration you can set it back to 25, which should be enough to handle the small fluctuations that happen during a simulation.

[zhenyuwei99]

Great! I think all my questions have been solved. Really thanks!