External Code Replacing Total PE Calculation

Jacob_Faibussowitsch · May 28, 2020, 1:02am

Hello All,

I am writing a code coupling a few libraries together (among them LAMMPS) to perform MD. I would like to replace how LAMMPS computes the total (and by extension local) potential energy for system relaxation with my own implementation. I have the following restrictions in no particular order:

Driver code is in C.
LAMMPS is not the driver.
User “Setup” for the code is minimal. Assuming the other libraries are separately built, they should simply be able to point code to the libraries and then compile and run. I want as much as possible to avoid a lengthy and error-prone configure stage. Ideally all I would need is a LAMMPS executable.

My current workflow is as follows:

I call MPI_init from the driver, perform setup and process command line options.
The driver then calls lammps_open passing sanitized command line options to LAMMPS.
The user is required to also pass a .in LAMMPS-style script via command line, so that use of my code is identical to using plain LAMMPS.
After parsing the .in script, I pass all lines before either a “run” or “minimize” command to LAMMPS as normal.
Based on these commands I single-step LAMMPS simulation (by calling the same run or minimize command only with time steps of 1) in a loop, performing my total energy calculation at every timestep.

At this point I would like to pass this back value to LAMMPS (or better yet, have LAMMPS call my code at every timestep itself to save me having to search for suitable places to hook in skipping steps 4 and 5), but I have not found a way to do this. Initially I tried replacing the compute->compute_scalar() function pointer for the internal compute thermo_pe that is created for every instance, but I have since learned that monkey patching is illegal in C and C++. I have also looked at fix_external, which is precisely the plug-and-play structure I need, but as far as I can tell I can only influence forces from a fix and provide a component to the energy summation rather than outright replace it.

Short of shipping my code with its own compute_pe.cpp and replacing the existing file, then recompiling LAMMPS (which incidentally would violate restriction 3) I am a bit out of ideas. Worst case I suppose I would have to loosen the last restriction somewhat, but perhaps there is a clever solution for this.

akohlmey · May 28, 2020, 3:10am

jacob,

unless you are changing the energy as part of something that is called in the “middle” of a timestep, it will have no effect whatsoever.

so using fix_external to execute a callback function is the way to go.

from this point onward, however, there is currently no way to continue with C, you need to dig fairly deep into the inner workings of LAMMPS and use some C++ code. i am currently working on a significant refactoring and improved documentation for the library interface(s). I think it should be possible to update the library interface to not only set the callback function but also call the functions to set the per atom or global energy or stress values for the fix, if combined with first defining a suitable compute, you can use the call back to “set” the energy by subtracting the existing energy contributions from the value you want to set.

back to the “hackish” C++ variant. the fewest lines of C++ code, would be to directly overwrite the energy in the corresponding accumulator variables.
assuming you have only a pair style, those are LAMMPS_NS::Pair::eng_vdwl and LAMMPS_NS::Pair::eng_coul. if you have other styles contributing to the (potential) energy, those need to be overwritten as well.

please see below of a short test code. if you put the mycallback() function into a separate file, then this is the only file that needs to be compiled with a C++ compiler. it has to use C bindings to be compatible with the callback signature.

#include “library.h”
#include “lammps.h”
#include “force.h”
#include “pair.h”
#include “fix_external.h”

#include <mpi.h>
#include <stdio.h>
#include <string.h>

extern “C” void mycallback(void *ptr, long int timestep,
int nlocal, int *ids,
double **x, double **fexternal);

void mycallback(void *ptr, long int timestep,
int nlocal, int *ids,
double **x, double **fexternal)
{
LAMMPS_NS::LAMMPS *lmp = (LAMMPS_NS::LAMMPS )ptr;
double my_pe = (500.0 - timestep)/100.0; / new energy: (500-step)/100 */

lmp->force->pair->eng_vdwl = my_pe;
lmp->force->pair->eng_coul = 0.0;
}

int main(int argc, char **argv)
{
void *handle;
int i;

MPI_Init(&argc, &argv);
lammps_open(0, NULL, MPI_COMM_WORLD, &handle);
lammps_file(handle,“in.ar.real”);
lammps_command(handle,“thermo_style custom step pe”);
lammps_command(handle,“thermo 1”);
lammps_command(handle,“fix ext all external pf/callback 1 1”);
lammps_set_fix_external_callback(handle,(char *)“ext”,&mycallback,handle);
lammps_command(handle,“run 0 post no”);

lammps_command(handle,“run 1 pre no post no”);
lammps_command(handle,“run 1 pre no post no”);
lammps_command(handle,“run 1 pre no post no”);
lammps_command(handle,“run 1 pre no post no”);
lammps_command(handle,“run 1 pre no post no”);

lammps_close(handle);
MPI_Finalize();
return 0;
}

this produces the following output:

LAMMPS (5 May 2020)
using 1 OpenMP thread(s) per MPI task
Lattice spacing in x,y,z = 5.79518 5.79518 5.79518
Created orthogonal box = (0 0 0) to (28.9759 28.9759 28.9759)
1 by 1 by 1 MPI processor grid
Created 500 atoms
create_atoms CPU = 0.00173871 secs
Neighbor list info …
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 9.8112
ghost atom cutoff = 9.8112
binsize = 4.9056, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Setting up Verlet run …
Unit style : real
Current step : 0
Time step : 11.1903
Per MPI rank memory allocation (min/avg/max) = 2.644 | 2.644 | 2.644 Mbytes
Step Temp PotEng Press v_temp v_epair v_emol v_etotal v_press
0 117.7 -1.6612661 -1909.5509 1 -7.1026383 0 -5.6056383 -5.1203128
10 87.369977 -1.5728967 -1169.6414 0.74231077 -6.7248204 0 -5.6135812 -3.1363029
20 42.567295 -1.4427006 -150.87379 0.36165926 -6.1681752 0 -5.6267713 -0.40455638
30 55.130978 -1.480902 -387.17817 0.46840253 -6.3315028 0 -5.6303042 -1.0381883
40 55.054202 -1.4807485 -401.72653 0.46775023 -6.3308469 0 -5.6306248 -1.0771986
50 56.873955 -1.4860029 -428.9126 0.48321117 -6.3533113 0 -5.6299442 -1.1500959
60 58.33701 -1.490161 -458.23636 0.49564154 -6.3710892 0 -5.6291138 -1.2287253
70 61.29671 -1.4991528 -539.72484 0.52078768 -6.4095331 0 -5.629914 -1.4472304
80 63.214984 -1.504992 -594.34987 0.53708567 -6.4344983 0 -5.630481 -1.5937032
90 61.936907 -1.5013008 -569.13985 0.5262269 -6.4187169 0 -5.6309552 -1.5261045
100 62.20662 -1.5023046 -585.49121 0.52851844 -6.4230083 0 -5.6318162 -1.5699494
Loop time of 0.0972865 on 1 procs for 100 steps with 500 atoms

Performance: 993.809 ns/day, 0.024 hours/ns, 1027.892 timesteps/s
99.0% CPU use with 1 MPI tasks x 1 OpenMP threads

akohlmey · May 28, 2020, 2:22pm

FWIW, I have now a version that can be called from the C library interface using C code only.

it requires some additions to the library interface to be able to call the (public C++) functions inside of fix_external() to set the energy.
to be able to change the energy from the fix, it needs to define a compute pe instance that computes all potential energy except that from fixes (so any fix that also changes the energy cannot be used with this). the callback will then look something like this:

void mycallback(void *handle, long int timestep,
int nlocal, int *ids,
double **x, double *fexternal)
{
double my_pe = (500.0 - timestep)/100.0; / new energy: (500-step)/100 */
double *pe = (double *)lammps_extract_compute(handle,“pe”,LMP_STYLE_GLOBAL,LMP_TYPE_SCALAR);
if (pe) my_pe -= *pe;
lammps_fix_external_set_energy_global(handle, “ext”, my_pe);
}

This change is currently in a feature branch of mine and it will not be ready for the next patch release of LAMMPS (planned for early next week), but hopefully for the one after that (in 3-4 weeks).

HTH,
Axel.

Jacob_Faibussowitsch · May 28, 2020, 5:55pm

Hello Axel,

Thank you for the in-depth explanations! I will implement both examples you have given, and swap them around once your feature branch is merged.

it requires some additions to the library interface to be able to call the (public C++) functions inside of fix_external() to set the energy.

Why not extend the same fix_external functionality and create for example a compute_external? I’m not hugely familiar (more than just some minor digging) with LAMMPS internals so perhaps there are some irreconcilable differences between the two…

This change is currently in a feature branch of mine and it will not be ready for the next patch release of LAMMPS (planned for early next week), but hopefully for the one after that (in 3-4 weeks).

Is this on a pushed branch? If so I can check it out and test my code against it locally.

akohlmey · May 28, 2020, 6:15pm

Hello Axel,

Thank you for the in-depth explanations! I will implement both examples you have given, and swap them around once your feature branch is merged.

it requires some additions to the library interface to be able to call the (public C++) functions inside of fix_external() to set the energy.

Why not extend the same fix_external functionality and create for example a compute_external?

there is no point to it. compute pe is not actually doing a computation, but rather returns a sum of data that is already collected internally in different places, all of those are (derived) class instances that may be present or not. the main point about computes is that they need a “consumer” to flag to the real computing part of LAMMPS when to collect that data. if it is not required (e.g. for thermo output), then the calculations of those terms are skipped. this applies to energies and the virial (which is used to compute pressure).

I’m not hugely familiar (more than just some minor digging) with LAMMPS internals so perhaps there are some irreconcilable differences between the two…

there simply is no obvious value to that proposition. what you are describing you are trying do is not exactly a common situation.

replacing energy that is consistent with the computed forces with energy that is not, is quite unusual. the available functionality is quite suitable for the common case where some part of the force calculation or some kind of constraint or other force manipulation is delegated to some external code, which then returns modifications to the force (and energy, and stress as needed).

This change is currently in a feature branch of mine and it will not be ready for the next patch release of LAMMPS (planned for early next week), but hopefully for the one after that (in 3-4 weeks).

Is this on a pushed branch? If so I can check it out and test my code against it locally.

this is work in progress and only partially tested and may change (and break) any time as i am continuing to refactor the library interfaces and rework its documentation. if you still want access, please send me an e-mail off list and i can point you into the right direction, but i will ignore any questions about the code in that branch until it is merged.

axel.