How to compute a custom force for each particle MD simulation

[Adrien-Berard]

Hello, I am new to hoomd-blue and I would like to compute forces and energy for each particles that are not in the provided functions. I did the research online but could not find any working answers (maybe I did it wrong). Can anyone help me please ? Here is the code that I can provide and I would like to run a simple MD simulation of a polymer :

Code

import numpy as np

import hoomd

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)
simulation.create_state_from_gsd(filename='random.gsd')

from scipy.special import lambertw

def compute_b_lambertw(l0):
b_solution = lambertw(-2 * np.exp(-2)) / 2
return b_solution.real

Exponential_force_particle_energy function

def exponential_force_particle_energy(positions, bonds, l0, b, particle_index):
num_particles = len(positions)
force = np.zeros(3, dtype=np.float32)
energy = 0.0

force -= 4 / l0 * np.sum(
    np.sign(positions[:, np.newaxis] - positions[particle_index]) *
    np.exp(-4 * np.abs(positions[:, np.newaxis] - positions[particle_index]) / l0) / l0,
    
)
energy += np.sum(
    np.exp(-4 * np.abs(positions[:, np.newaxis] - positions[particle_index]) / l0)
)

# Calculate bonding forces and energy
for bond in bonds:
    if particle_index in bond:
        partner = bond[1] if particle_index == bond[0] else bond[0]
        bond_force = 4 * b / l0 * np.sign(positions[particle_index] - positions[partner]) * \
                     np.exp(-4 * b * np.abs(positions[particle_index] - positions[partner]) / l0) / l0
        force += bond_force
        energy -= np.sum(np.exp(-4 * b * np.abs(positions[particle_index] - positions[partner]) / l0))

return force, energy

Example usage

particle_index = 0 # Replace with the desired particle index
force_on_particle, energy_on_particle = exponential_force_particle_energy(positions, bonds, l0, b, particle_index)
print(f"Force on particle {particle_index}: {force_on_particle}")
print(f"Energy on particle {particle_index}: {energy_on_particle}")

Define your potential parameters

l0 = 1.0
b = compute_b_lambertw(l0)

Define system parameters

dimension = 3
num_particles = 10
diffusion_coefficient = 1.0
time_step = 1e-4
num_steps = 1000

Example bonded pairs

bonds = [(i, (i + 1) % num_particles) for i in range(num_particles - 1)]

Example initial positions

positions = np.random.rand(num_particles, dimension)

def exponential_force_energy_all(positions, bonds, l0, b):
num_particles = len(positions)
forces = np.zeros((num_particles, 3), dtype=np.float32)
energies = np.zeros(num_particles, dtype=np.float32)

# Calculate forces and energies for each particle
for i in range(num_particles):
    force, energy = exponential_force_particle_energy(positions, bonds, l0, b, i)
    forces[i] = force
    energies[i] = energy

return forces, energies

Example usage

forces_all_particles, energies_all_particles = exponential_force_energy_all(positions, bonds, l0, b)

Print forces and energies for all particles

for i in range(num_particles):
print(f"Force on particle {i}: {forces_all_particles[i]}")
print(f"Energy on particle {i}: {energies_all_particles[i]}")

integrator = hoomd.md.Integrator(dt = 0.005, integrate_rotational_dof=False, forces=exponential_force_energy_all(positions,bonds,l0,b)[0], constraints=None, methods=None, rigid=None)

returns : ValueError: Validation failed: Error raised in conversion: Value [14.122901 17.188808 16.567728] is not convertable into any of these types (<class 'hoomd.md.force.Force'>,)

Thanks for your help,

Adrien Berard

[joaander]

As the error message states, members of forces must be subclasses of hoomd.md.force.Force and override the appropriate methods. The guide How to apply arbitrary forces in MD in the documentation gives an overview of methods to apply arbitrary forces. Follow the links for details on each method.

[Adrien-Berard]

Ok thank you for that answer.
I tried the 3 options written on your link but I don't succeed to make any of them work. Here is my code for the last option where I create a custom force in python :

Python Code

import numpy as np

import hoomd
import gsd.hoomd
import matplotlib

%matplotlib inline
matplotlib.style.use('ggplot')
import matplotlib_inline

matplotlib_inline.backend_inline.set_matplotlib_formats('svg')

MD Simulation using Hoomd Blue

from scipy.special import lambertw

def compute_b_lambertw(l0):
b_solution = lambertw(-2 * np.exp(-2)) / 2
return b_solution.real
import numpy as np
import matplotlib.pyplot as plt

Parameters for the exponential potential

l0 = 1.0
b = compute_b_lambertw(l0)

Function to compute the potential energy

def exponential_potential(r, l0, b):
return np.exp(-4 * r / l0) - np.exp(-4 * b * r / l0)

Function to compute the force from the potential

def exponential_force(r, l0, b):
return 4 / l0 * (np.exp(-4 * b * r / l0) - np.exp(-4 * r / l0))

Generate a range of distances

r = np.linspace(0, 10 * l0, 500)

Calculate the potential energy and force for each distance

U_exponential = exponential_potential(r, l0, b)
F_exponential = exponential_force(r, l0, b)

Plot the exponential potential and force on the same graph

fig, ax1 = plt.subplots(figsize=(8, 6))

color = 'tab:red'
ax1.set_xlabel('Distance (r)')
ax1.set_ylabel('Potential Energy (U)', color=color)
ax1.plot(r, U_exponential, color=color, label='Exponential Potential')
ax1.tick_params(axis='y', labelcolor=color)

ax2 = ax1.twinx()
color = 'tab:blue'
ax2.set_ylabel('Force (F)', color=color)
ax2.plot(r, F_exponential, color=color, label='Force')
ax2.tick_params(axis='y', labelcolor=color)

fig.tight_layout()
plt.title('Exponential Potential and Force')
plt.show()

import numpy as np

def exponential_repulsion_force_energy(positions, l0):
num_particles = len(positions)
repulsion_forces = np.zeros((num_particles, 3), dtype=np.float32)
repulsion_energies = np.zeros(num_particles, dtype=np.float32)

for i in range(num_particles):
    repulsion_forces[i] -= 4 / l0 * np.sum(
        np.sign(positions[:, np.newaxis] - positions[i]) *
        np.exp(-4 * np.abs(positions[:, np.newaxis] - positions[i]) / l0) / l0,
    )
    repulsion_energies[i] = np.sum(np.exp(-4 * np.abs(positions[:, np.newaxis] - positions[i]) / l0))

return repulsion_forces, repulsion_energies

import numpy as np

def exponential_bond_force_energy(positions, bonds, l0, b):
num_bonds = len(bonds)
bond_forces = np.zeros((num_bonds, 3), dtype=np.float32)
bond_energies = np.zeros(num_bonds, dtype=np.float32)

for i, bond in enumerate(bonds):
    particle_i, particle_j = bond
    bond_direction = positions[particle_j] - positions[particle_i]
    bond_distance = np.linalg.norm(bond_direction)

    force_magnitude = 4 * b / l0 * np.exp(-4 * b * bond_distance / l0) / l0
    bond_force = force_magnitude * bond_direction / bond_distance

    bond_forces[i] = bond_force
    bond_energies[i] = -np.sum(np.exp(-4 * b * bond_distance / l0))

return bond_forces, bond_energies

Example usage

positions = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]], dtype=np.float32)
bonds = [(0, 1), (1, 2)]
l0 = 1.0
b = 0.5

repulsion_forces, repulsion_energies = exponential_repulsion_force_energy(positions, l0)
bond_forces, bond_energies = exponential_bond_force_energy(positions, bonds, l0, b)

for i in range(len(bonds)):
print(f"Bond {bonds[i]}:")
print(f"Bonding Forces: {bond_forces[i]}")
print(f"Bonding Energy: {bond_energies[i]}")
print()

Define your potential parameters

l0 = 1.0
b = 0.5

Define system parameters

dimension = 3
num_particles = 5
diffusion_coefficient = 1.0
time_step = 1e-4
num_steps = 1000

Example bonded pairs

bonds = [(i, (i + 1) % num_particles) for i in range(num_particles - 1)]

Example initial positions

positions = np.random.rand(num_particles, dimension)

repulsion_forces, repulsion_energies = exponential_repulsion_force_energy(positions, l0)
bond_forces, bond_energies = exponential_bond_force_energy(positions, bonds, l0, b)
bonds, bond_forces
class MyCustomForce(hoomd.md.force.Custom):
def init(self):
super().init(aniso=False)

def set_forces(self,positions,bonds,l0,b):
    with self.cpu_local_force_arrays as arrays:
        arrays.force[:] = exponential_bond_force_energy(positions,bonds,l0,b)[0] + exponential_repulsion_force_energy(positions,l0)[0]
        arrays.torque[:] = 0
        arrays.potential_energy[:] = exponential_bond_force_energy(positions,bonds,l0,b)[1] + exponential_repulsion_force_energy(positions,l0)[1]
        arrays.virial[:] = np.arange(6)[None, :]

frame = gsd.hoomd.Frame()

Place a polymer in the box.

frame.particles.N = 5
frame.particles.position = [[-2, 0, 0], [-1, 0, 0], [0, 0, 0], [1, 0, 0],
[2, 0, 0]]
frame.particles.types = ['A']
frame.particles.typeid = [0] * 5
frame.configuration.box = [20, 20, 20, 0, 0, 0]

Connect particles with bonds.

frame.bonds.N = 4
frame.bonds.types = ['A-A']
frame.bonds.typeid = [0] * 4
frame.bonds.group = [[0, 1], [1, 2], [2, 3], [3, 4]]

with gsd.hoomd.open(name='molecular.gsd', mode='x') as f:
f.append(frame)
forces = MyCustomForce()

Perform the MD simulation.

sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)
sim.create_state_from_gsd(filename='molecular.gsd')
langevin = hoomd.md.methods.Langevin(filter=hoomd.filter.All(), kT=1.0)
integrator = hoomd.md.Integrator(dt=0.005,
methods=[langevin]
)
gsd_writer = hoomd.write.GSD(filename='molecular_trajectory.gsd',
trigger=hoomd.trigger.Periodic(1000),
mode='xb')
sim.operations.integrator = integrator
sim.run(0)
integrator.forces.append(forces.set_forces(positions = positions,bonds = bonds,l0 = l0,b = b))
sim.operations.integrator = integrator
sim.operations.writers.append(gsd_writer)
sim.run(10e3)

I got this error :

Error

---

DataAccessError Traceback (most recent call last)
Cell In[73], line 13
11 sim.operations.integrator = integrator
12 sim.run(0)
---> 13 integrator.forces.append(forces.set_forces(positions = positions,bonds = bonds,l0 = l0,b = b))
14 sim.operations.integrator = integrator
15 sim.operations.writers.append(gsd_writer)

Cell In[70], line 6
5 def set_forces(self,positions,bonds,l0,b):
----> 6 with self.cpu_local_force_arrays as arrays:
7 arrays.force[:] = exponential_bond_force_energy(positions,bonds,l0,b)[0] + exponential_repulsion_force_energy(positions,l0)[0]
8 arrays.torque[:] = 0

File ~/opt/anaconda3/envs/myenv2/lib/python3.12/site-packages/hoomd/md/force.py:187, in Force.cpu_local_force_arrays(self)
183 raise RuntimeError("Cannot enter cpu_local_force_arrays context "
184 "manager inside another local_force_arrays "
185 "context manager")
186 if not self._attached:
--> 187 raise hoomd.error.DataAccessError("cpu_local_force_arrays")
188 return hoomd.md.data.ForceLocalAccess(self, self._simulation.state)

DataAccessError: The property cpu_local_force_arrays is unavailable until the simulation runs for 0 or more steps.

[joaander]

This is too hard to read. Can you please post a single question with your code and output formatted in code blocks.