Logging custom forces

[Beuteltee]

Hello everyone,
since the code iam using is quite long and seperated over multiple files iam reproducing the important lines.
So my problem is that i have a simulation running with multiple custom forces using the hoomd.md.force.Custom class and now i am willing to log the forces each for its own. To create a GSD file where all the different forces per time step are listed. I implemented a logger follwing the documentation(https://hoomd-blue.readthedocs.io/en/latest/tutorial/02-Logging/02-Saving-Array-Quantities.html) in the following way:

class custom_force(hoomd.md.force.Custom):

(here just how the force is defined according to the documentation. This already works since the simulation is running and yielding the right output.)

def force_writer(config:dict, cell: hoomd.md.nlist, filename='forces.gsd'):

logger=hoomd.logging.Logger(categories=['scalar','particle'])
logger.add(custom_force(config,cell),quantities=['forces'])

gsd_force_writer= hoomd.write.GSD(
    filename=filename,
    trigger=hoomd.trigger.And([
        hoomd.trigger.Periodic(config['recording_interval']),
        hoomd.trigger.Before(config['recording_interval'] * (
                config['n_burn_in_steps'] + config['n_recorded_steps_agent_observables'] + 1))
    ]),
    mode='wb',
    filter=hoomd.filter.All(),
    logger=logger
)
return gsd_force_writer

then there is a simulation class where i add and call the writer like this
self.gsd_force_writer = setup_force_writer(self.config,
setup_nlist(self.config),
f"{self.output_path_no_ext}_force.gsd")

self.simulation.operations.writers.append(self.gsd_force_writer)

self.gsd_force_writer.flush()

the GSD file actually is generated but when i try
traj = gsd.hoomd.open("forces.gsd", mode="r")

traj[0].log.keys()

i get

dict_keys([])

as output. Somehow this also happens for the trajectory GSD file, however when i convert it to a H5 file the file includes the trajectory.
Iam not sure now if i have to add a loggable to my custom force or if the force defined in the class like

class custom_force(hoomd.md.force.Custom):
def set_forces(self, timestep):
with self.cpu_local_force_arrays as arrays:
arrays.force[:] =...

is already logged as quantitiy=["forces"]. Would be nice if somone knows anything about it. It would also be helpful if i need some other type of command to verify if there is something written in the GSD file.

[joaander]

You don't provide a full minimal working example, so I can only guess. You appear to be taking all the needed steps:

1. Including 'particle' in the logger categories.
2. Adding the force quantity to the logger.
3. Including the logger in the GSD writer.

There are many moving parts to your code. To troubleshoot, you need to break the problem in to smaller parts.

1. Check the return value of logger.log() and ensure that it includes the force.
2. Check the return value of your custom force's forces property (I don't see what variable name you assign that to, you appear to pass an unnamed temporary to logger.add).
3. Try adding simulation to the logger and see if its loggables appear in the GSD file.

These tests will help you narrow down whether the problem is with the Logger, CustomForce or the GSDWriter. If you still can't solve the problem, then you will know exactly what to include (and NOT to include) in a minimal working example.

A minimal working example should use hoomd.util.make_example_simulation, demonstrate which ONE of the three checks above fails. It should be a single copy and pasted code block (formatted as code) with as few lines as possible and print() to show the achieved result. I should be able to copy and paste it into a file and run it with python script.py. You do not need to repeat your entire custom force code (should the CustomForce be needed in a MWE at all) - a class with pass in set_forces is sufficient to demonstrate whether there is a problem logging or not.

[Beuteltee]

Thanks for your quick response. I was working a bit with it and, since the code where i am trying to add the functionality to log and save the forces of the simulation is not from me, i decided to test out logging according to the documentation. I came up with the following code:

import gsd.hoomd 
import hoomd
import matplotlib
import numpy`
#matplotlib inline
#matplotlib.style.use('ggplot')
import matplotlib_inline

matplotlib_inline.backend_inline.set_matplotlib_formats('svg')

class test_force(hoomd.md.force.Custom):
    def __init__(self):
        super().__init__(aniso=True)

    def set_forces(self, timestep):
        with self.cpu_local_force_arrays as arrays:
            arrays.force[:] = -5

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

integrator = hoomd.md.Integrator(dt=0.005)
cell = hoomd.md.nlist.Cell(buffer=0.4)
lj = hoomd.md.pair.LJ(nlist=cell)
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)
lj.r_cut[('A', 'A')] = 2.5

integrator.forces.append(lj)
integrator.forces.append(test_force())

nvt = hoomd.md.methods.ConstantVolume(
    filter=hoomd.filter.All(), thermostat=hoomd.md.methods.thermostats.Bussi(kT=1.5)
)
integrator.methods.append(nvt)
simulation.operations.integrator = integrator
simulation.run(0)

print(test_force().loggables['forces'])

logger = hoomd.logging.Logger(categories=['particle'])
logger.add(lj, quantities=['energies', 'forces'])
logger.add(test_force(), quantities=['energies','forces'])

gsd_writer = hoomd.write.GSD(
    filename='trajectory.gsd',
    trigger=hoomd.trigger.Periodic(10000),
    mode='xb',
    filter=hoomd.filter.All(),
)
simulation.operations.writers.append(gsd_writer)

gsd_writer.logger = logger

simulation.run(100000)

gsd_writer.flush()

traj = gsd.hoomd.open('trajectory.gsd', mode='r')

print(traj[0].log.keys())

print(traj[-1].log['particles/md/pair/LJ/forces'][0:10]) 

which creates the following output:

*Warning*: Forces provide torques, but integrate_rotational_dof isfalse.
particle
*Warning*: Forces provide torques, but integrate_rotational_dof isfalse.
dict_keys(['particles/md/pair/LJ/energies', 'particles/md/pair/LJ/forces'])
[[ 1.75859327e-13  1.84297022e-13 -2.58681965e-13]
 [ 4.65183447e-14 -3.24573701e-13  9.21485110e-14]
 [-6.75015599e-14  1.29007915e-13 -5.75095527e-14]
 [-2.57571742e-14  2.85105273e-13  2.28927988e-13]
 [ 1.62578284e-13  3.64604180e-14  6.72795153e-14]
 [-1.06581410e-14 -6.57252031e-14  5.28466160e-14]
 [ 1.80300219e-13 -6.88782364e-13  3.50830476e-14]
 [-7.67164110e-13 -8.65751915e-13  2.26985097e-13]
 [-2.68673972e-14 -1.44551038e-13 -5.30686606e-14]

So it shows that the test_force class is working, since print(test_force().loggables['forces']) creates the particle output and also the force array at the end is different to when i leave test_force out of the simulation. But the Logging still doesnt work as i expected it. I was hoping to see the ['energies'] and ['forces'] in the dict_keys output but they are missing. So to refer to your answer i think test 1 and 3 are working but test 2 fails so it might be a problem on how i add the force to the logger.

[joaander]

Store an instance of your test force class in a variable and reuse it:

test_force_instance = test_force()
integrator.forces.append(test_force_instance)
# ...
print(test_force_instance.loggables['forces'])
# ...
logger.add(test_force_instance, quantities=['energies','forces'])

The values are not logged because the specific test force instance that you pass to the logger is not connected to the simulation and has no state to operate on. You should get an error if you attempt to compute forces on such an instance:

test_force().forces # error: simulation must be run 1 or more steps

[Beuteltee]

Thank you a lot that fixed my problem!