How to apply distance constraints on a chain of spring-beads?

[yliu6677]

I have a chain of spring-beads. I want to make the spring rigid while keep the angular freedom. I saw the class hoomd.md.constrain.Distance(tolerance=0.001) might be able to achieve this. But I didn't see an example like the rigid constraint. I was wondering how to use this distance constraint to a specific group of particles in the simulation.

Alternatively, I can set the spring constant such that the springs are so stiff if this distance constrain cannot be applied.

[joaander]

constrain.Distance applies the constraints listed in the constraints portion of the simulation state which can be set in a GSD file

The unit test provides a working example:

hoomd-blue/hoomd/md/pytest/test_constrain_distance.py

Lines 64 to 75 in 41af37e

	for x in x_coords:
	for i, y in enumerate(y_coords):
	position.append([x, y, 0])
	if i & 1:
	constraint_values.append(bead_spacing)
	tag = len(position) - 1
	constraint_groups.append([tag, tag - 1])
	s.particles.position[:] = position
	s.constraints.N = len(constraint_values)
	s.constraints.value[:] = constraint_values
	s.constraints.group[:] = constraint_groups

hoomd-blue/hoomd/md/pytest/test_constrain_distance.py

Lines 127 to 133 in 41af37e

	d = hoomd.md.constrain.Distance()
	sim = simulation_factory(polymer_snapshot_factory())
	integrator = hoomd.md.Integrator(dt=0.005)
	nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
	integrator.methods.append(nve)
	integrator.constraints.append(d)

[yliu6677]

I made two minimum examples:

import itertools
import math

import gsd.hoomd # 3.4.2
import hoomd # 5.0.0
import matplotlib.pyplot as plt
import numpy as np

parameters = {'dt': 0.0001,
              'kT': 1.,
              'rna_monomer_radius': 1,
              'rna_monomer_mass': 1.,
              'rna_length': 16,
              'box_length': 100,
              'k_s_rna_nomoner_harmonic': 100,
              'brownian_gamma': 1.,
              'brownian_gamma_r': (1., 1., 1.),
              'gsd_writer_trigger_period': 1,
              'gsd_writer_mode': 'wb',
              'simulation_run_time': 1000,
              'rna_chain_monomer_constrain_tol':0.01
              }

rna_chain_init_radius = parameters['rna_monomer_radius'] / np.sin(np.pi / parameters['rna_length'])
rna_chain_monomer_constrain_spacing = 2 * parameters['rna_monomer_radius']

initial_frame = gsd.hoomd.Frame()

initial_frame.particles.N = parameters['rna_length']
initial_frame.particles.position = [[rna_chain_init_radius * np.cos(p * 2 * np.pi / parameters['rna_length']), rna_chain_init_radius * np.sin(p * 2 * np.pi / parameters['rna_length']), 0] for p in np.arange(parameters['rna_length'])]
initial_frame.particles.orientation = [(1, 0, 0, 0)] * parameters['rna_length']
initial_frame.particles.types = ["rna_monomers"]
initial_frame.particles.typeid = [0] * parameters['rna_length']
initial_frame.configuration.box = [parameters['box_length'], parameters['box_length'], parameters['box_length'], 0, 0, 0]

initial_frame.constraints.N = parameters['rna_length'] - 1
initial_frame.constraints.value = (parameters['rna_length'] - 1) * [rna_chain_monomer_constrain_spacing]
initial_frame.constraints.group = [[i, (i + 1)] for i in np.arange(parameters['rna_length'] - 1)]

with gsd.hoomd.open(name="rna_init.gsd", mode='w') as f:
    f.append(initial_frame)
    f.flush()
    f.close()

simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)
simulation.create_state_from_gsd(filename='rna_init.gsd')
brownian = hoomd.md.methods.Brownian(filter=hoomd.filter.All(), kT=parameters['kT'], default_gamma=parameters['brownian_gamma'], default_gamma_r=parameters['brownian_gamma_r'])
rna_chain_monomer_constrain = hoomd.md.constrain.Distance(tolerance=parameters['rna_chain_monomer_constrain_tol'])
integrator = hoomd.md.Integrator(dt=parameters['dt'], integrate_rotational_dof=False, forces=[], methods=[brownian], constraints=[rna_chain_monomer_constrain])
simulation.operations.integrator = integrator
gsd_writer = hoomd.write.GSD(trigger=hoomd.trigger.Periodic(parameters['gsd_writer_trigger_period']), filename='dynamics_rna.gsd', mode=parameters['gsd_writer_mode'])
simulation.operations.writers.append(gsd_writer)
simulation.run(10)

gsd_writer.flush()

and

import itertools
import math

import gsd.hoomd # 3.4.2
import hoomd # 5.0.0
import matplotlib.pyplot as plt
import numpy as np

parameters = {'dt': 0.0001,
              'kT': 1.,
              'rna_monomer_radius': 1,
              'rna_monomer_mass': 1.,
              'rna_length': 16,
              'box_length': 100,
              'k_s_rna_nomoner_harmonic': 100,
              'brownian_gamma': 1.,
              'brownian_gamma_r': (1., 1., 1.),
              'gsd_writer_trigger_period': 1,
              'gsd_writer_mode': 'wb',
              'simulation_run_time': 1000,
              'rna_chain_monomer_constrain_tol':0.01
              }

rna_chain_init_radius = parameters['rna_monomer_radius'] / np.sin(np.pi / parameters['rna_length'])
rna_chain_monomer_constrain_spacing = 2 * parameters['rna_monomer_radius']

initial_frame = gsd.hoomd.Frame()

initial_frame.particles.N = parameters['rna_length']
initial_frame.particles.position = [[rna_chain_init_radius * np.cos(p * 2 * np.pi / parameters['rna_length']), rna_chain_init_radius * np.sin(p * 2 * np.pi / parameters['rna_length']), 0] for p in np.arange(parameters['rna_length'])]
initial_frame.particles.orientation = [(1, 0, 0, 0)] * parameters['rna_length']
initial_frame.particles.types = ["rna_monomers"]
initial_frame.particles.typeid = [0] * parameters['rna_length']
initial_frame.configuration.box = [parameters['box_length'], parameters['box_length'], parameters['box_length'], 0, 0, 0]

initial_frame.constraints.N = parameters['rna_length'] - 1
initial_frame.constraints.value = (parameters['rna_length'] - 1) * [rna_chain_monomer_constrain_spacing]
initial_frame.constraints.group = [[i, (i + 1)] for i in np.arange(parameters['rna_length'] - 1)]

with gsd.hoomd.open(name="rna_init.gsd", mode='w') as f:
    f.append(initial_frame)
    f.flush()
    f.close()

simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)
simulation.create_state_from_gsd(filename='rna_init.gsd')
integrator = hoomd.md.Integrator(dt=0.005)
nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
integrator.methods.append(nve)
rna_chain_monomer_constrain = hoomd.md.constrain.Distance(tolerance=parameters['rna_chain_monomer_constrain_tol'])
integrator.constraints.append(rna_chain_monomer_constrain)
simulation.operations.integrator = integrator
simulation.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=0.1)
gsd_writer = hoomd.write.GSD(trigger=hoomd.trigger.Periodic(parameters['gsd_writer_trigger_period']), filename='dynamics_rna.gsd', mode=parameters['gsd_writer_mode'])
simulation.operations.writers.append(gsd_writer)
cell = hoomd.md.nlist.Cell(buffer=0.4)
lj = hoomd.md.pair.LJ(nlist=cell)
lj.params[("rna_monomers", "rna_monomers")] = dict(epsilon=1, sigma=1)
lj.r_cut[("rna_monomers", "rna_monomers")] = 2 ** (1 / 6)
integrator.forces.append(lj)
simulation.run(10000)

gsd_writer.flush()

The first one has fatal problems which I didn't resolve by having a smaller dt, gamma or kT. The second one is fine for the currently choice of constrain tolerance. Can you tell the problem? I wanted to have both the Brownian motion and the constraint.

[joaander]

The Lagrange multiplier scheme implemented in constrain.Distance cannot be made compatible with the Brownian integration method. It is valid only for Velocity-Verlet, which in principle means that it should only work with ConstantVolume and no thermostat. ConstantVolume or ConstantPressure with a thermostat interoperate with constrain.Distance well enough in practice as the errors remain bounded.

If you have other particles in your system, you can use ConstantVolume with no thermostat on your constrained chains and place them in contact with a thermostatted solvent. The solvent-chain interactions will effectively give you Brownian motion.

If you require Brownian integration on your chain particles, than harmonic springs are your only option. Be warned that the Brownian integration method is only accurate to $O(\delta t)$, so you will need a very small step size to resolve large spring constants. The Langevin integrator (with error $O(\delta t^2)$ is much more efficient if your model can tolerate nonzero inertia.

[yliu6677]

Is the difference between Brownian and Langevin is that Brownian does not explicitly depends on the mass of the particle, does not consider the drag force and uses a different integration method. I tried Langevin and now the chain appears to evolve under the constraint, with a much slower motion though. Is it possible to make this Langevin motion a little more violent? Is it consistent with simulation.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1). It appears to be unaffected by this line.

How does the thermal torque apply on the particle? Suppose I have a single spherical particle. Will the torque change the status of this particle? How about a tetrahedron particle or a center particle in a general rigid body?

In this part of the model, I just wanted to add some random motion to the chain, instead of leaving it static. It is desired to make the chain rigid constraint, whereas a spring-bead chain is also acceptable. There will be another type of solute that interacts with this chain. That rigid body will have Brownian motion, so they can meet and interact with each other. The solvent will not be shown explicitly but through the modification of the interaction forces between the solutes.