Assistance Developing Compute

LAMMPS LAMMPS Mailing List Mirror
1

I’m currently using LAMMPS as part of a proton transfer simulation, and in the course of the simulation, I need to extract the energy of the system at two different configurations (based on proton location) during each time step. The way I was doing this originally, which appears to work, was to use the LAMMPS library interface with the following sequence of commands:

run 1 (to propagate system)
reset proton coordinates to state 1
run 0 (to get potential energy)
reset proton coordinates to state 2
run 0
revert to original coordinates and continue simulation

Of course, this required each step to be propagated separately, and the method is relatively slow. In the interest of speed and simplifying the interface, I’ve been working to develop a compute that would perform the same calculation during each time step automatically, without requiring additional communication from an external program. Schematically, the relevant bit of my existing compute code looks something like:

Compute * potential = modify->compute[pe_index];

move_atom(x1);

first_energy = potential->compute_scalar();

move_atom(x2);

second_energy = potential->compute_scalar();

Although I’ve checked that the functions are moving the proper atom to the correct location, it appears that the potential energy is not updated between calls, so I wind up with both energy values identical for every step, which is obviously wrong, and conflicts with known results using the previous library method. I’ve tried doing a few things to force a recalculation of the energy, but haven’t had any luck so far and I was wondering if someone could help me figure out how to tell the pe compute to reset between my evaluation at states 1 and 2?

I’ve looked through the relevant source files for the pe compute, and the verlet and integrate modules, but I’m not entirely clear on how these manage to require an update on each time step.

Thanks,
Thomas Allen

2

Potential energies are collected during the force calculation, if requested. Thus it has to be signaled beforehand whether enegies are needed. So your principal idea cannot work.

You can try to speed up your library version if you can skip the preparation by using the run pre no option.

Otherwise you would have to use partitions like NEB to parallelize the multiple force computes.

Axel.

ps: there should be a discussion about how to implement a simple EVB model in the mailing list archives, which seems to go into a similar direction.

3

correct location, it appears that the potential energy is not updated
between calls, so I wind up with both energy values identical for every

You need to trigger a force calculation. This can be done with
force->pair->compute(1, 0);

You may also have to force a kspace update using:
force->pair->compute(1, 0);

A faster solution would be calculating the differential energy. I.e.
avoid recalculating parts of the box that do not change as a function
of the proton location by iterating only over the proton neighborhood.
This solution will however have to be implemented for every potential
class you want this to work for.

4

sorry the second call would have to be:
force->kspace->compute(1, 0);