Calculating the energy of a trajectory

[Ramon Crehuet]

Dear HOOMD users,

I am quite new to HOOMD. I have two (apparently) very simple questions, which I am not able to solve by reading the documentation.

1. Is it possible to store the potential energy of the configuration that is saved in a GSD trajectory?  dynamic does not recognize the quantity 'potential_energy', even though analyze.log does. I know I could save the energy in a text file with the logger, I was just wondering if it could be stored in the GSD file.
   
2. If I load a trajectory with gsd.hoomd.open('trajectory.gsd', 'rb'), can I calculate the energy of the each frame? I guess I also need to define a system (with the same number of particles and the interactions), but once this is done, how can I assign the coordinates in the trajectory to the system?

Context explanation:

I plan to run long trajectories keeping only uncorrelated frames in the GSD. I would like to calculate the energy of these frames when I slightly modify the parameters of the interaction potential. That is why I wanted to store the initial energy of the frames and why I need to re-calculate the energy of the frames with the new parameters.

Thanks for your time.

Best regards,

Ramon

[Joshua A. Anderson]

Ramon,

Thanks for asking. The just released HOOMD-blue v3.0.0-beta.1 has a completely new API that we designed specifically with these types of use-cases in mind. There aren't specific tutorials for these yet in the documentation. Here are complete scripts to demonstrate using a LJ liquid.

You can save the potential energy of the system in a GSD file using hoomd.md.compute.ThermodynamicQuantities and hoomd.logging.Logger:
```
import hoomd

# configure the Simulation
gpu = hoomd.device.GPU()
sim = hoomd.Simulation(device=gpu)
sim.create_state_from_gsd(filename='random.gsd')

# set up the LJ liquid MD operations
integrator = hoomd.md.Integrator(dt=0.005)
cell = hoomd.md.nlist.Cell()
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)
nvt = hoomd.md.methods.NVT(kT=1.5, filter=hoomd.filter.All(), tau=1.0)
integrator.methods.append(nvt)
sim.operations.integrator = integrator

# log the thermodynamic quantities to an output GSD file
log = hoomd.logging.Logger()
thermo = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())
sim.operations.computes.append(thermo)
log += thermo   # this logs all loggable quantities: potential_energy, kinetic_energy, etc ...
# log.add(thermo, quantities=['potential_energy'])  # this logs only potential_energy

gsd = hoomd.write.GSD(filename='trajectory.gsd',
                      trigger=hoomd.trigger.Periodic(10000),
                      overwrite=True)
gsd.log = log
sim.operations.writers.append(gsd)

sim.run(1e5)
```

OVITO can read these quantities from the GSD file. You can access them with the GSD python package:
```
import gsd.hoomd

with gsd.hoomd.open('trajectory.gsd') as f:
    for frame in f:
        print(frame.log['md/compute/ThermodynamicQuantities/potential_energy'])
```

To answer your second question, this script uses HOOMD v3 to compute the potential energy of a system over frames in a GSD file:
```
import hoomd

gpu = hoomd.device.GPU()

for frame in range(10):
    # configure the Simulation
    sim = hoomd.Simulation(device=gpu)
    sim.create_state_from_gsd(filename='trajectory.gsd', frame=frame)

    # set up the LJ liquid MD operations
    integrator = hoomd.md.Integrator(dt=0.005)
    cell = hoomd.md.nlist.Cell()
    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)
    sim.operations.integrator = integrator

    # compute the system energy
    thermo = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())
    sim.operations.computes.append(thermo)

    # prepare the computation, but don't actually run any time steps
    sim.run(0)

    # get the potential energy
    print(thermo.potential_energy)
```

In the last example, the Integrator is necessary because it is the object that sums the forces to get a total potential energy. Future beta releases may allow you to compute (and/or log) the energy term from a Force object that is not associated with the integrator equations of motion.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/410fcf26-c475-4528-a96c-81d3c0c9102an%40googlegroups.com.

[Ramon Crehuet]

Hi Joshua,

That was so very helpful! Thank you.

I was also reading the paper in Scipy 2020 (http://conference.scipy.org/proceedings/scipy2020/pdfs/brandon_butler.pdf) and from there I thought of a different way to achieve my second question: assigning just the particle position. I thought it could be faster, because I do not build the simulation at each step. I also thought it was more clear in my case, as I change the parameters once before reading the trajectory. My attempt is the following:

new_energy = []
# configure the Simulation
sim = hoomd.Simulation(device=cpu)
sim.create_state_from_gsd(filename='trajectory.gsd', frame=0)

# set up the LJ liquid MD operations
integrator = hoomd.md.Integrator(dt=0.005)
cell = hoomd.md.nlist.Cell()
lj = hoomd.md.pair.LJ(nlist=cell)

lj.params[('A', 'A')] = dict(epsilon=1.0, sigma=1)

lj.r_cut[('A', 'A')] = 2.5
integrator.forces.append(lj)
sim.operations.integrator = integrator

snap = sim.state.snapshot

# compute the system energy
thermo = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())
sim.operations.computes.append(thermo)

with gsd.hoomd.open('trajectory.gsd') as f:
    for frame in f:

        snap.particles.position[:,:] = frame.particles.position
        sim.state.snapshot = snap

        # prepare the computation, but don't actually run any time steps

        sim.run(1)

        # get the potential energy

        new_energy.append(thermo.potential_energy)

It seems to work, but I have some doubts. First, I need to run 1 step because with sim.run(0) the calculated energies are all those of the first frame. Why?

And second, if I set sim.run(1) I get the warning:

*Warning*: integrate.mode_standard: No integration methods are set, continuing anyways.

In fact, it seems that no steps are run, as the energy is the same as the one from the GSD file. Still, I would like to remove this warning and have a nicer code.

All the best,

Ramon

El dia divendres, 16 d’octubre de 2020 a les 15:01:36 UTC+2, Joshua A. Anderson va escriure:

[Ramon Crehuet]

OK, I think I understand the warning (it was obvious): there isn't any integrator method. If I add:

langevin = hoomd.md.methods.Langevin(hoomd.filter.All(), kT=1., seed=42)
integrator.methods.append(langevin)

The warning disappears. The problem now is that it actually performs 1 intergrator step, and the energy vary. That is not what I want. Is there a workaround?

Thanks again,

Ramon

El dia dimarts, 20 d’octubre de 2020 a les 13:04:00 UTC+2, Ramon Crehuet va escriure:

[Joshua A. Anderson]

Ramon,

Since all you are doing is computing the energy, you can safely ignore the warning. The warning is intended to help users that are running a simulation notice when they forget to add any methods. You can work around the warning by adding an NVE method with dt=0.

The reason the energy doesn't change with run(0) is because HOOMD assumes that quantities like forces and energies don't change until the timestep number changes. The first run(0) will compute the initial net force and energy and the next will not be computed until step 1. You can work around this by querying `lj.energy` directly as it should recompute after you assign the snapshot.

Another way to write this would be to implement your frame loading code in a custom updater and call sim.run() once for the number of frames in the file.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

 1. Is it possible to store the potential energy of the configuration that is saved in a GSD trajectory?  dynamic does not recognize the quantity 'potential_energy', even though analyze.log does. I know I could save the energy in a text file with the logger, I was just wondering if it could be stored in the GSD file.

 2. If I load a trajectory with gsd.hoomd.open('trajectory.gsd', 'rb'), can I calculate the energy of the each frame? I guess I also need to define a system (with the same number of particles and the interactions), but once this is done, how can I assign the coordinates in the trajectory to the system?

Context explanation:

I plan to run long trajectories keeping only uncorrelated frames in the GSD. I would like to calculate the energy of these frames when I slightly modify the parameters of the interaction potential. That is why I wanted to store the initial energy of the frames and why I need to re-calculate the energy of the frames with the new parameters.

Thanks for your time.

Best regards,

Ramon

--

You received this message because you are subscribed to the Google Groups "hoomd-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/410fcf26-c475-4528-a96c-81d3c0c9102an%40googlegroups.com.

 --

You received this message because you are subscribed to the Google Groups "hoomd-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/8e5a6ce9-bd90-4c34-93ca-3934f06ef07an%40googlegroups.com.

[Ramon Crehuet]

Hi Joshua,

Thanks again. Setting dt=0 removed the integrator warning but raised another warning: A timestep of less than 0.0 was specified. Which is strictly not

true as the timestep is exactly equal to zero and not less than zero. Maybe you could correct the code now that it is still in beta version.

But the use of lj.energy worked smoothly.

Cheers,

Ramon

El dia dimarts, 20 d’octubre de 2020 a les 14:06:08 UTC+2, Joshua A. Anderson va escriure: