Logging energies for harmonic external potentials

[edwinb-ai]

Hello,

I am trying to compute some free energies of solids, and I have made the following code that should log the positions of the particles during a production run, as well as the energies from the harmonic external potential that link the particles to its lattice sites. However, when I check the output of the trajectory, there are no energies logged.

I have checked that the external potential does have the energy field to log, but the resulting trajectory is empty whenever I check the keys() of the resulting file.

Here is the code. I really am not sure what I am doing wrong, any help is be appreciated.

from pathlib import Path

import freud
import gsd.hoomd
import hoomd
import numpy as np

PROJECT_DIR = Path().resolve()


def simulation(spring_constant, rep):
    # Create the initial lattice, in this case, FCC
    fcc = freud.data.UnitCell.fcc()
    box, positions = fcc.generate_system(num_replicas=4)
    N_particles = positions.shape[0]
    print(f"Number of particles: {N_particles}")
    # Define the density to work with and compute the size of the box
    low_density = 1.04086
    low_L = np.cbrt(N_particles / low_density)
    # Re-scale the positions of the particles to the new box
    positions *= [low_L, low_L, low_L] / box.L
    # Define here the reference positions and orientations for the
    # harmonic potential implementation
    reference_positions = np.copy(positions)
    reference_orientations = [(1, 0, 0, 0)] * N_particles
    # Then we deal with the initialization of the system for the simulation
    frame = gsd.hoomd.Frame()
    frame.particles.N = N_particles
    frame.particles.position = positions
    frame.particles.orientation = np.copy(reference_orientations)
    # There are two types of particles: mobile and immobile, the immobile
    # particle is the reference particle (carrier according to Vega & Noya)
    frame.particles.types = ["hs", "hs_no_motion"]
    frame.particles.typeid = np.zeros((N_particles,))
    # This is the ID for the carrier particle for Vega-Noya method
    frame.particles.typeid[0] = 1
    frame.configuration.box = [low_L, low_L, low_L, 0, 0, 0]

    # Create a directory with the information of this state point
    data_path = PROJECT_DIR.joinpath(
        f"Npart={N_particles}",
        f"data_{rep}",
        f"k={spring_constant:.6g}_npart={N_particles}",
    )
    data_path.mkdir(parents=True, exist_ok=True)
    with gsd.hoomd.open(name=data_path.joinpath("lattice.gsd"), mode="w") as f:
        f.append(frame)

    # We start by defining the simulation object
    device = hoomd.device.CPU()
    rng = np.random.default_rng()
    sim = hoomd.Simulation(device=device, seed=rng.integers(2**16))
    sim.create_state_from_snapshot(frame)
    # We define the integrator
    mc = hoomd.hpmc.integrate.Sphere(default_d=0.01)
    mc.shape["hs"] = dict(diameter=1.0)
    mc.shape["hs_no_motion"] = mc.shape["hs"]
    # Set movement to zero for the single particle
    mc.a["hs_no_motion"] = 0.0
    mc.d["hs_no_motion"] = 0.0
    sim.operations.integrator = mc

    # During equilibration, we need the harmonic constraints activated
    harmonic = hoomd.hpmc.external.Harmonic(
        reference_positions,
        reference_orientations,
        k_translational=spring_constant,
        k_rotational=1.0,
        symmetries=[[1, 0, 0, 0]],
    )
    mc.external_potential = [harmonic]
    device.notice(harmonic.loggables)

    # * Equilibrate the system
    # Add a displacement tuner with a target acceptance ratio of 0.4
    tune = hoomd.hpmc.tune.MoveSize.scale_solver(
        moves=["d"], target=0.4, trigger=hoomd.trigger.Periodic(50)
    )
    sim.operations.tuners.append(tune)
    device.notice("Starting equilibration run...")
    sim.run(1e5)
    device.notice(f"Acceptance rate: {mc.translate_moves[0] / sum(mc.translate_moves)}")
    sim.operations.tuners.remove(tune)

    # * Production stage
    # Add logging
    logger = hoomd.logging.Logger(categories=["scalar"])
    logger.add(harmonic, quantities=["energy"])

    gsd_writer = hoomd.write.GSD(
        filename=data_path.joinpath("trajectory.gsd"),
        trigger=hoomd.trigger.Periodic(100),
        mode="wb",
        logger=logger,
        filter=hoomd.filter.All(),
    )
    sim.operations.writers.append(gsd_writer)
    device.notice("Collecting data...")
    sim.run(1e5)
    device.notice("Done!")


if __name__ == "__main__":
    # Test with a spring constant of 2.0
    simulation(2.0)

[edwinb-ai]

I found the issue, which is in this line

mc.external_potential = [harmonic]

I replaced it with

sim.operations.integrator.external_potentials = [harmonic]

so that the integrator already assigned was updated correctly.

[joaander]

In your case, mc and sim.operations.integrator are references to the same object.

The effective change was external_potential = [harmonic] -> external_potentials = [harmonic]. The former was (almost) valid in HOOMD-blue 4.x. The latter is valid in HOOMD-blue 5.x (https://hoomd-blue.readthedocs.io/en/v5.1.1/migrating.html#migrating-to-hoomd-blue-5).

Specifically, in HOOMD-blue 4.x, external_potential is a single potential (functioning code would be external_potential = harmonic). In 5.x, you can set multiple external potentials. Now you provide a list of potentials in the external_potentials attribute. It is a "feature" of Python that callers can add new class attributes, so mc.external_potential = <anything> should not generate a warning or error even though the HOOMD-blue code ignores it.

[edwinb-ai]

Thanks for the clear explanation! Makes a lot of sense.