Setting patch width for patchy potentials

[ccasert]

Hi,

I'm having problems with setting the patch width correctly for patchy potentials in MD simulations. Regardless of the value of alpha I choose in the envelope_params, I get the same results - even when I set it to zero. Below is a minimal script that shows the energy is non-zero even when alpha = 0. What am I doing wrong here? I'm using v5.1.0

Thanks!

import hoomd
import gsd.hoomd
import numpy as np
import itertools
num_particles = 100
theta = np.random.random(num_particles)*2*np.pi-np.pi
orientations = [tuple([np.cos(the/2),0,0,np.sin(the/2)]) for the in theta]

L = 10
x = np.linspace(-L / 2, L / 2, 10, endpoint=False)
positions = np.array(list(itertools.product(x, repeat=2)))
positions = np.hstack([positions, np.zeros((positions.shape[0], 1))])

frame = gsd.hoomd.Frame()
frame.particles.N = num_particles
frame.particles.typeid = [0] * num_particles
frame.particles.position = positions
frame.particles.orientation = orientations
frame.configuration.box = [L, L, 0, 0, 0, 0]

frame.particles.types = ['A']
frame.particles.moment_inertia = [0,0, 1/8] * num_particles

with gsd.hoomd.open(name="init.gsd", mode='w') as f:
    f.append(frame)

dt = 1e-3

device = hoomd.device.auto_select()
sim = hoomd.Simulation(device=device, seed = np.random.randint(0, 65535))
sim.create_state_from_gsd("init.gsd")

pot_attr = hoomd.md.pair.aniso.PatchyExpandedLJ(nlist = hoomd.md.nlist.Cell(buffer=0.4))
envelope_params=dict(alpha=0, omega=50)
lj_params = dict(epsilon=10, sigma = 1/2, delta = 1 - .5*2**(1/6))

pot_attr.params[('A', 'A')] = dict(pair_params=lj_params, envelope_params=envelope_params)
pot_attr.directors['A']     =  [(np.cos(np.pi/6),np.sin(np.pi/6),0), (np.cos(-np.pi/6),np.sin(-np.pi/6),0),(np.cos(np.pi/6 + np.pi),np.sin(np.pi/6+np.pi),0), (np.cos(-np.pi/6+np.pi),np.sin(-np.pi/6+np.pi),0) ]
pot_attr.r_cut[('A','A')] = 1.5



nvt = hoomd.md.methods.ConstantPressure(
filter=hoomd.filter.All(),
thermostat=hoomd.md.methods.thermostats.Bussi(kT=1),
tauS=100*dt,
S=1,
gamma = 0,
couple="xy",
)

sim.operations.integrator = hoomd.md.Integrator(dt=dt, methods = [nvt], forces = [pot_attr],integrate_rotational_dof=True )
sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1)

logger = hoomd.logging.Logger(categories=['scalar'])
logger.add(pot_attr, quantities=['energy'])

file = open("log.txt", mode='w', newline='\n')

table = hoomd.write.Table(output=file,trigger=hoomd.trigger.Periodic(period=100), logger=logger)
sim.operations.writers.append(table)

sim.run(int(1e4))

[cbkerr]

Hi @ccasert, thanks for the report and steps to reproduce!

I see two parts in your comment: energies not changing and energies being non-zero even when alpha is 0. I couldn't reproduce the first part. For alpha=np.pi/3, I get energy around -5.3e+03 and for alpha=np.pi/2 energy of -1.0e+04.

Alpha=0 is not designed to turn off a patch, but I agree that this is unintuitive behavior. To turn off a patch for a pair interaction, set its r_cut to 0. Does that work for you?

Having a non-zero energy is a byproduct of the design of the normalization. I agree that this is confusing. We only validate that the user-provided alpha is within 0 to pi inclusive. It might make sense to change that to an exclusive interval.

[joaander]

In addition to what @cbkerr said, I notice that you don't have an isotropic repulsive core. Without it, your particles can pass through each other when the patches are not aligned. If they later rotate to align patches, then they might also interact on the r < delta side of the potential.