Requesting a full neighbor list in C++ code using LAMMPS as a library

LAMMPS LAMMPS Development
,
1

Hi all!

I am developing a C++ code that uses LAMMPS as a library.

Once I initialize LAMMPS and run for a few timesteps, I can easily access the pair-wise neighbor list requested by the interatomic potential through the relevant pointer as follows:

NeighList *mylist = lmp->force->pair->list;

Here, I face a couple of issues. Firstly, the neighbor list requested by the interatomic potential I am using (lj/cut/coul/cut) is a “half list”. However, I definitely need a “full list” for the code I am developing.

Secondly, I am using a “special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0” command in my LAMMPS input file. This means that the neighbor list requested by the pair style will not include “bonded atoms”. Of course, there is a way of accessing the bond list using something like that:

int **bondlist = lmp->neighbor->bondlist;

Generally, it would be really helpful if I could somehow request a new neighbor list from within the C++ script that is “full” and includes bonded neighboring atoms. This way I could simply use a single data container, let LAMMPS deal with the scattering of data across different processes, and nicely iterate through my neighbor list to retrieve the necessary data.

Is there any way I could request a new “full” neighbor list that also includes bonded interactions, from within the C++ code? Can I “delete/remove” the neighbor list once I am done using it and perform an additional run?

Here is an example code:

//System Headers
#include <iostream>
#include <vector>
#include <fstream>
#include <cstdio>
#include <chrono>
#include <cstdlib>
#include <cstring>
#include <mpi.h>

//LAMMPS Headers
#include "./libraries/LAMMPS_source/lmppath.h"
#define QUOTE_(x) #x
#define QUOTE(x) QUOTE_(x)
#include QUOTE(LMPPATH/src/lammps.h)
#include QUOTE(LMPPATH/src/input.h)
#include QUOTE(LMPPATH/src/atom.h)
#include QUOTE(LMPPATH/src/library.h)
#include QUOTE(LMPPATH/src/neighbor.h)
#include QUOTE(LMPPATH/src/neigh_list.h)
#include QUOTE(LMPPATH/src/force.h)
#include QUOTE(LMPPATH/src/pair.h)
#include QUOTE(LMPPATH/src/compute.h)
#include QUOTE(LMPPATH/src/pointers.h)

using namespace LAMMPS_NS;
using namespace std;

void setupMPI( int *narg , char ***arg , int *nprocs , int *me , double &time_start ){
	
	MPI_Init( narg , arg );
	MPI_Comm_rank( MPI_COMM_WORLD , me );
	MPI_Comm_size( MPI_COMM_WORLD , nprocs );
	time_start = MPI_Wtime( );
	
}

void initializeLMP( LAMMPS **lmp ){

	*lmp = new LAMMPS( 0 , NULL , MPI_COMM_WORLD );
	
}


void readLMPinput( const std::string &lmp_input , LAMMPS *lmp ){
	
	lmp->input->file( lmp_input.c_str( ) ); //LAMMPS constructor only uses c style strings so 
	
}

void finalizeMPI( int *me , const double &time_start ){

	double time_end = MPI_Wtime( );
	
	printf( "Elapsed time is %f\n (process %i) \n", time_end - time_start , *me ); 
	MPI_Finalize();
	
}

void deleteArraysAndObj( LAMMPS **lmp ){

	delete *lmp;
	
}

int main( int narg , char **arg )
{
	//Declare variables
	int me , nprocs;
	double time_start;
	string lmp_input = arg[ 2 ];
	LAMMPS *lmp = NULL;
	
	//Prepare
	if ( narg != 3 ) {
		printf( "Syntax Error. Input command should be in the following form: mpirun -np N ./main -in in.lammps\n" );
		exit( 1 );
	}
	setupMPI( &narg , &arg , &nprocs , &me , time_start );
	initializeLMP( &lmp );
	readLMPinput( lmp_input , lmp );
	
	
	//Main
	lmp->input->one("run 1000"); //RUN LAMMPS FOR A FEW STEPS
	
	
	//Request a "full" neighbor list that includes bonded interactions -> HOW???
	//Loop over neighbor list in parallel and extract data/perform operations
	//Delete neighbor list
	
	lmp->input->one("run 1000"); //RUN LAMMPS OF A FEW MORE STEPS
	
	
	
	
	//Finalize
	deleteArraysAndObj( &lmp );
	finalizeMPI( &me , time_start );
}

Thank you in advance for your help :slight_smile:

Kind regards,
Stavros

2

Since you are using the C library interface elsewhere, you should be using the corresponding library interface call. 1.1.7. Neighbor list access — LAMMPS documentation
That is more consistent and safer. Also, direct access is not very future-proof. Protecting internal data structures of LAMMPS classes from global access has been on the agenda for making LAMMPS more reliable and fit for the future for quite some time now (and already implemented for some parts).
So using documented accessor functions is the way to go instead of accessing data members of classes directly.

There is no way to build a neighbor list from an external code, since there always must be a requestor class that is derived from known base style in LAMMPS (Pair, Compute, Fix, Bond, Command, etc.). But you could change your pair style to:

pair_style hybrid/overlay lj/cut/coul/cut <cutoff> zero <cutoff> full
pair_coeff * * zero

and then adapt the pair_coeff commands for the rest to include lj/cut/coul/cut as needed for hybrid styles. Pair style zero will do no computation, but does request a neighbor list. The “full” setting will request a full neighbor list. And so you can look up the pair style neighbor list for style “zero” to get access to a full neighbor list.

Technically, it would also be possible to create a custom pair style of your own that then is instantiated outside of the Force class and contains a custom neighbor list request, but that is rather messy and requires a more in-depth understanding of LAMMPS’ internals. The hybrid/overlay approach is thus much safer, specifically when the neighbor list is looked up through the C library interface.

Another option could be to implement your feature as a custom Command style in LAMMPS. Those can request and explicitly trigger the build of an occasional neighbor list (unlike pair styles which request perpetual neighbor lists).

Neighbors are only excluded, if their special bonds scale factor is exactly zero. If you would use 1.0e-100 instead, the result would be the same (the scale factor is far beyond the precision limits of double precision floating point math), but the neighbors are included. But you have to note that the 1-2 neighbors are flagged with extra bits on the neighbor’s atom index and those must be masked off or else you get a segmentation fault. For details, please see the discussion of how to implement the Pair::compute() method or a pair style here: 4.7.1. Writing new pair styles — LAMMPS documentation

Neighbor lists are automatically destroyed at the beginning of a run and then the neighbor list code figures out how to build the neighbor lists according to the active requests in the most efficient way.

2 Likes
3

Hi @akohlmey.

Thank you so much for your kind help and your fast response. Everything is really helpful, really appreciated!

I attempted to change the pair style to:

pair_style hybrid/overlay lj/cut/coul/cut 12.0 zero 12.0 full

However, I get this error when I try to run LAMMPS:

ERROR: Illegal pair_style command (src/pair_zero.cpp:95)
Last command: pair_style hybrid/overlay lj/cut/coul/cut 12.0 zero 12.0 full

I don’t get this error when I remove the “full” argument at the end of the command.

Is there any quick way to use a hybrid/overlay with the zero full command?

Again, thank you for your kind help.

4

If you look at the online documentation you will find a note New in version 3Nov2022 next to the full keyword. So you are likely using an older version of LAMMPS.

When doing development with LAMMPS we recommend to stay up-to-date, so that porting efforts remain minimal if there is ever something coming out of that that should be contributed back to LAMMPS. Also, you have the latest bugfixes and improvements. If you are adventurous, this means following the “develop” branch, or if you are more cautious, then this means following the “release” branch". Following “stable” or “maintenance” or just using a copy of the downloaded source instead of a checkout of a github fork, will immediately put you at a disadvantage.

1 Like
5

Hi @akohlmey!

I compiled my code with the “release” branch and everything works nicely now.

I used the following code (as in LAMMPS documentation) to mask off the extra bits of the neighbor’s atom index to avoid the segmentation fault:

j &= NEIGHMASK;

Thank you for your help!

Stavros

1 Like