Question about velocity-verlet integrator and python API

LAMMPS LAMMPS Development
1

Hi, everyone,

i have a question about how to build my integrator using python.

my team’ve develop a new integration algorithm and want to make it by lammps python api. therefore, to build a framework for algorithm migration, we firstly use python package of lammps to implement velocity-verlet integrator (VV) (lammps integrator by default).

the VV integrator made by ourselves works under NVE ensemble but does not work under NVT ensemble. here are the core function.

def gradientFunction(Lammps: IPyLammps | PyLammps, Position: np.ndarray) -> np.ndarray:
    """
    define the gradient function
    """
    llnp = Lammps.lmp.numpy

    x: np.ndarray = llnp.extract_atom('x')
    x[:] = Position

    before_id = llnp.extract_atom('id').copy()

    Lammps.run(0, 'pre yes post no');

    after_id = llnp.extract_atom('id')

    assert (before_id == after_id).all(), 'array has changed !!!!!'

    force: np.ndarray = llnp.extract_atom('f')

    return force.copy()

def VelocityVerlet(Lammps: PyLammps | IPyLammps, basis_timestep: float, ntimestep: int, disable_tqdm: bool = False):

    def excute(Lammps: PyLammps | IPyLammps, Dt: float):
        # get the initial state of system
        lmpX, lmpV = Lammps.lmp.numpy.extract_atom('x'), Lammps.lmp.numpy.extract_atom('v')
        lmpMass, atomtype = Lammps.lmp.numpy.extract_atom('mass'), Lammps.lmp.numpy.extract_atom('type')

        # copy variables
        x_0, v_0 = Lammps.lmp.numpy.extract_atom('x').copy(), Lammps.lmp.numpy.extract_atom('v').copy()
        mass = rearrange(lmpMass[atomtype], "l -> l 1")
        
        # keep the atoms of boundary stationary
        matrix = np.ones_like(x_0)
        matrix[atomtype < 3] = 0.
        v_0[atomtype < 3] = 0.

        ftm2v = ftm2v_coeff[Lammps.system.units]
        massinv = matrix * ftm2v / mass
        Dt_half = Dt * 0.5

        # get gradient at t by using Lammps API and compute v(t + Dt/2)
        force = gradientFunction(Lammps, x_0)
        v_ht = v_0 + massinv * force * Dt_half

        # compute r(t + Dt)
        x_t = x_0 + v_ht * Dt

        # get gradient at t + Dt by using Lammps API and compute v(t + Dt)
        force = gradientFunction(Lammps, x_t)
        v_t = v_ht + massinv * force * Dt_half

        # put x_t, v_t into source address
        lmpX[:] = x_t; lmpV[:] = v_t

    with tqdm(total = ntimestep, desc = 'vv Iteration: ', leave = False, position = 1, disable = disable_tqdm) as vv_bar:
        for _ in range(ntimestep):
            excute(Lammps, basis_timestep)
            vv_bar.update()

    return

here are the lammps code

def create_simulation(timestep = 0.2, cmdargs = None, num_threads: int = 1, ensemble: str = 'nve') -> IPyLammps:

    lmp = lammps(cmdargs=cmdargs)
    
    MDsimulation = PyLammps(ptr = lmp)

    MDsimulation.enable_cmd_history = True
    if num_threads > 1:
        MDsimulation.package(f"omp {num_threads} neigh yes")
        MDsimulation.suffix('omp')
    
    MDsimulation.units('real')
    MDsimulation.atom_style('full')

    MDsimulation.pair_style('lj/cut/coul/long 12.0')
    MDsimulation.bond_style('harmonic')
    MDsimulation.angle_style('harmonic')
    MDsimulation.dihedral_style('opls')
    MDsimulation.improper_style('cvff')

    MDsimulation.dielectric(1.0)
    MDsimulation.pair_modify('mix arithmetic')
    MDsimulation.special_bonds('lj/coul 0.0 0.0 1.0')
    
    MDsimulation.read_data('input/data.lammps')

    MDsimulation.atom_modify('sort 0 0.0') # turn off sort algorithm

    MDsimulation.set('type 1 charge -0.55')
    MDsimulation.set('type 2 charge 1.1')

    MDsimulation.group('zeo type 1 2 ')
    MDsimulation.group('indole type 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18')
    
    MDsimulation.kspace_style('pppm', 1e-4)

    MDsimulation.neighbor('10.0 bin')
    MDsimulation.neigh_modify('every 1 delay 0 check yes exclude molecule/intra zeo')

    MDsimulation.delete_bonds('zeo multi')

    MDsimulation.velocity('indole create 700.0 902144 dist gaussian')
    MDsimulation.velocity('zeo set 0.0 0.0 0.0')
    MDsimulation.velocity('indole scale 700.0')

    if ensemble == 'nvt':
        MDsimulation.fix('1 indole nvt temp 700.0 700.0 1.0')
        # MDsimulation.fix("1 indole temp/rescale 1 700.0 700.0 0.01 1.0")
        # MDsimulation.fix("2 indole nve")
    elif ensemble == 'nve':
        MDsimulation.fix('1 indole nve')

    MDsimulation.compute('1 indole temp')

    MDsimulation.thermo_modify('lost/bond ignore')
    MDsimulation.thermo_style("custom step time spcpu temp press pe ke enthalpy evdwl emol ecoul etotal")

    MDsimulation.timestep(timestep) # attn set timestep

    # initialize system state
    MDsimulation.run(0, 'pre yes post no')

    return MDsimulation

here 's the energy variation graph over time in the NVE ensemble. it works.

velocity-verlet integrator in the NVE ensemble
velocity-verlet integrator in the NVE ensemble1920×525 76.3 KB

and here’s energy variation graph over time in the NVT ensemble.

velocity-verlet integrator in the NVT ensemble
velocity-verlet integrator in the NVT ensemble1920×525 90.1 KB

how can i implement the function for correct data?
Thank you so much

2

Doing any serious programming with PyLammps is a bad idea.

Ignoring lost atoms or bonds is almost always a bad idea. They mean that there is something wrong.

A neighbor list skin of 10.0 is significant waste of CPU time.

Back when I was a graduate student my adviser would call such panels “stamp collection” meaning that lots of graphs were collected, but no commentary and discussion of what was important and no selection of which were the graphs indicating a noteworthy behavior and where.

Overall, it is not clear what exactly you are doing and how. There doesn’t seem to be a connection between the first script and the second script.

If you want to implement / test a new integrator for LAMMPS in Python, then your approach is misguided. It does not account for particles moving between subdomains.

Instead you want to compile LAMMPS with the PYTHON package enabled and then implement it via the fix python/move command — LAMMPS documentation

This has been used with success to study time integration algorithms with Python code. Check out this publication: https://doi.org/10.1080/00268976.2019.1649493

3

thank you for your tips.

And i mean, is there a method that i can directly call the LAMMPSFixMove (in fix python/move command) class in main python script to move atoms without (re)compling it? because recompling is a bit agonizing for debug process.

another question is that i have no idea for the meaning of parameters such as vflag, ilevel and iloop in LAMMPSFixMove class and how to use these parameters.

4

No.

Adding the PYTHON package is a one time step. Once the fix command is available, you only need to modify python code (and you can also load the python module to access internal data of LAMMPS).

Even if you are writing code in Python you need to understand some of the C++ code surrounding that you want to embed your calls into. Please have a look at, e.g.: 4.7. How a timestep works — LAMMPS documentation and 4.8.7. Writing a new fix style — LAMMPS documentation

That said, if you just want to use LAMMPS as a “force engine” and do everything else in Python, then both using the Python module and the PYTHON package are a good ideas. Instead you should look at the MDI software and package: https://docs.lammps.org/Packages_details.html#pkg-mdi and 8.1.7. Using LAMMPS with the MDI library for code coupling — LAMMPS documentation

If I remember correctly, MDI has a python module, so you can implement your timestepping entirely in Python and then concurrently set up your system to compute forces in LAMMPS and launch your python code and LAMMPS through the MDI software and have a two way client / server communication between your code and LAMMPS through the MDI protocol.
Of course, that also likely requires a recompile of LAMMPS to include the package (if it is not already included).

IMO, compilation of LAMMPS with CMake is very straightforward once you have understood a few CMake basics.

5

Hi, akohlmey

i’m thankful for giving these very useful tutorials to me. i decide to use PYTHON package for,which is convenient to transfer codes to other enviroments afterwards. But after reading it, i still have no idea how to make my new integrator.
In a timestep work. the following process of verlet integrator would be executed

loop over N timesteps:
  if timeout condition: break
  ev_set()
  
  fix->initial_integrate()
  fix->post_integrate()
  
  nflag = neighbor->decide()
  if nflag:
    fix->pre_exchange()
    domain->pbc()
    domain->reset_box()
    comm->setup()
    neighbor->setup_bins()
    comm->exchange()
    comm->borders()
    fix->pre_neighbor()
    neighbor->build()
    fix->post_neighbor()
  else:
    comm->forward_comm()
  
  force_clear()
  fix->pre_force()
  
  pair->compute()
  bond->compute()
  angle->compute()
  dihedral->compute()
  improper->compute()
  kspace->compute()
  
  fix->pre_reverse()
  comm->reverse_comm()
  
  fix->post_force()
  fix->final_integrate()
  fix->end_of_step()
  
  if any output on this step:
    output->write()

# after loop
fix->post_run()

but my algorithm may be as follows

loop over N timesteps:
  for i 1 to M: # add 
    if timeout condition: break
    ev_set()
  
    fix->initial_integrate()
    fix->post_integrate()
  
    nflag = neighbor->decide()
    if nflag:
      fix->pre_exchange()
      domain->pbc()
      domain->reset_box()
      comm->setup()
      neighbor->setup_bins()
      comm->exchange()
      comm->borders()
      fix->pre_neighbor()
      neighbor->build()
      fix->post_neighbor()
    else:
      comm->forward_comm()
  
    force_clear()
    fix->pre_force()
  
    pair->compute()
    bond->compute()
    angle->compute()
    dihedral->compute()
    improper->compute()
    kspace->compute()
  
    fix->pre_reverse()
    comm->reverse_comm()
  
    fix->post_force()
    fix->final_integrate()
    fix->end_of_step()
  
  if any output on this step:
    output->write()

# after loop
fix->post_run()

Unlike the respa method, the M may be greater than 30 (if i used respa method, i need to write a very long command including many 1) and our method does not involve sum_flevel_f() function (line 536 src/respa.cpp)

how should i write the process (or LAMMPSFixMove Class in fix python/move command) using PYTHON package and Python module?

thank you very much

6

You cannot. You would need to implement a new run style (in C++).

From your post it is not clear why you cannot just multiply the number of timesteps with M.
There don’t seem to be any additional changes to the flow of control.

[Update:]

To explain a bit more…
Please note that it is next to impossible to give meaningful advice on these questions with this extremely limited context you are providing. Your requirements as you have stated them are conflicting: a) you want to program entirely in Python, b) you want to use a workflow different from what is implemented in LAMMPS. You cannot have both. If you use fix python/move, you can write code only in Python, but then you must follow one of the available run styles (i.e. verlet of respa). If you don’t want this, you have to implement a new workflow. If you don’t want to do this in C++, the only option left to you is to reduce LAMMPS to a “force engine” and implement a client/server code using the MDI library as I had referred to before. In that case you have to write the entire flow of control as your MDI driver code where you then set up the model in LAMMPS and from then on keep sending LAMMPS new coordinates and retrieving forces from it and the rest has to be done by code written by you.

Bottom line: you either have to follow what LAMMPS offers and insert your python code where LAMMPS has designated interfaces or reduce LAMMPS’ role to that of a “force engine” or implement your changes in C++.