MuVT has problems keeping particles within the boundaries of a given triclinic box #2235
Description
Description
I am currently working on a project, where I want to simulate a Monte Carlo simulation of molecules on the surface of a triclinic crystal. In a first step, simplified the system to cuboids as particles and a triclinic box to find out via a GCMC (using MuVT) how many particles fit into the box. However, the particles receive coordinates that are far outside the box (at least 61x outside of it). This causes visualization tools (like Fresnel) to depict a highly distorted box. Testing the same simulation but with a cubic box, everything works fine. But with triclinic, it doesn't see where the box starts and ends.
Below is a the essential block from my code.
Script
import hoomd
import gsd.hoomd
import numpy as np
Lx = 143
Ly = 143 * np.sin(np.pi * (60 / 180)) # Tilt angle: 60°; for cubic test: 90°
Lz = 50
xy = 143 * np.cos(np.pi * (60 / 180))
frame = gsd.hoomd.Frame()
frame.configuration.box = [Lx, Ly, Lz, xy, 0, 0]
frame.particles.types = ["A"]
frame.particles.N = 1
frame.particles.position = [[0, 0, 0]]
frame.particles.typeid = [0]
with gsd.hoomd.open(name="src/convexPolyhedrons_in_triclinicBox.gsd", mode="w") as f:
f.append(frame)
# Start building the simulation.
simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed = 1)
simulation.create_state_from_gsd(filename="src/convexPolyhedrons_in_triclinicBox.gsd")
# Define ConvexPolyhedrons. Vertices list the coordinates of the corners of the Polyhedron.
mc = hoomd.hpmc.integrate.ConvexPolyhedron(default_d=0.1, default_a=0.1)
mc.shape["A"] = dict(
vertices=[
(+3.8, +3.45, +1.4),
(+3.8, +3.45, -1.4),
(+3.8, -3.45, +1.4),
(+3.8, -3.45, -1.4),
(-3.8, +3.45, +1.4),
(-3.8, +3.45, -1.4),
(-3.8, -3.45, +1.4),
(-3.8, -3.45, -1.4),
]
)
simulation.operations.integrator = mc
# Generate Hard Walls.
top = hoomd.wall.Plane(origin=(0,0,5), normal=(0,0,-1))
bottom = hoomd.wall.Plane(origin=(0,0,-5), normal=(0,0,1))
walls = [top, bottom]
wallPotential = hoomd.hpmc.external.wall.WallPotential(walls=walls)
simulation.operations.integrator.external_potentials = [wallPotential]
# Create the simulation updater.
muvt = hoomd.hpmc.update.MuVT(
trigger=hoomd.trigger.Periodic(10),
transfer_types=["A"],
)
muvt.fugacity["A"] = 100
simulation.operations.updaters.append(muvt)
# Create a tuner for the simulation.
daTuner = hoomd.hpmc.tune.MoveSize(
trigger = hoomd.trigger.Periodic(100),
target = 0.2,
moves = ['d', 'a'],
solver = hoomd.tune.ScaleSolver()
)
simulation.operations.tuners.append(daTuner)
# Set up the logger.
logger = hoomd.logging.Logger()
logger.add(mc, quantities=["translate_moves", "rotate_moves"])
# Set up a writer to transfer the data from the logger to a saved file.
writer = hoomd.write.GSD(
filename = "data/triclinic_trajectory.gsd",
trigger = hoomd.trigger.Periodic(1000),
mode = 'wb',
filter=hoomd.filter.All(),
logger = logger
)
simulation.operations.writers.append(writer)
# Equilibrate the simulation.
print("Starting simulation ...")
print("Phase 1: Filling Box. Tuning translational and rotational step size.")
simulation.run(100_000)
print("Phase 2: Perform the measurement.")
simulation.run(100_000)
final_N = simulation.state.N_particles
print(f"After performing the simulation, there are {final_N} particles in the system.")
print(f"The final frame looks as followed:")
render(simulation.state.get_snapshot(), "cPH_in_triBox_finished")Input files
For render, I used the method that is used in the hoomd-blue docs (slightly adjusted for my script, but essentially the same).
Output
My output can be seen in the images below. In the cubic system, the simulation is performed as intended, while in the simulation of the triclinic system, the box goes beyond the image boundaries there are barely any particles visibles. In both systems are about 110 - 120 particles at the end of the simulation.Expected output
Platform
Linux, CPU
Installation method
Conda-forge package
HOOMD-blue version
6.1.0
Python version
3.11.14