Extract (LJ) interatomic potential coefficients from LAMMPS object?

LAMMPS LAMMPS Development
,
1

Hi all!

I was wondering if there is any consistent way of extracting the (lj-cut-coul-cut) interatomic potential coefficients (sigma and epsilon) from a LAMMPS object.

I am using LAMMPS as a library for a C++ code. In the beginning, I used the following command to read the interatomic potential coefficients for any atom-type pair in the system:

lmp->input->file( "in.lmp" );

Then, I was able to acquire the interatomic coefficients as follows:

int types_num = *( int *)lammps_extract_global( lmp , const_cast<char*>( "ntypes" ) );

double **sigma_temp = ( double **)lmp->force->pair->extract( const_cast<char*>( "sigma" ) , 2 );

Afterward, I attempted to access the data using a double for loop:

	for( int i = 1; i < types_num + 1; ++i ){
			
		for( int j =1; j < types_num + 1; ++j ){

					printf( "Sigma%d%d=%f and sigma%d%d=%f \n" , i , j , sigma_temp[i][j] , j , i , sigma_temp[j][i] );

	
		}
			
			
	}

Here, I observed the following:

  1. For some pairs sigma_temp[i][j] has a non-zero value but sigma_temp[j][i] is zero
  2. For some other pairs sigma_temp[j][i] has a non-zero value but sigma_temp[i][j] is zero

I believe that overcoming issues 1 and 2 is not very challenging. One can simply create a new array and see if a pair has zero sigma_temp, and if so initialize it with the value of the inverse pair:(i.e., sigma[i][j] = sigma_temp[j][i] when sigma[i][j] = 0.0 and vice versa.). However, I was just generally wondering if retrieving the interatomic coefficients as above is thread-safe.

Also, I am aware that I can retrieve the interatomic potential coefficients as I read the file by writing a custom function (i.e., reading the input file line by line on proc 0, extracting the coefficients from there, and then broadcasting the information to all other procs). However, I was just wondering if there is any other thread-safe way of extracting epsilon or sigma directly from an object of the LAMMPS class (by using a lammps_extract_global function, for example).

Thank you in advance for your kind help :slight_smile:

Cheers,
Stavros

2

Which version of LAMMPS are you using?

This has nothing to do with this. Extracting the data is not subject to thread safety because it is a read-only operation for as long as you do not attempt it while LAMMPS is running in a different thread. If you intend to change the parameters, let me tell you right away: forget about it. It is not that simple. Just have a look at the code for fix adapt.

With the current version of LAMMPS, the matrix of per-type data should be fully symmetric after a first run. In general, you would also have to extract and inspect the “setflag” property.
only if setflag is set to 1, the parameter has been set from the input. The remaining settings will be generated during the init() phase of a run from the set mixing rule. Before, accessing those element results in undefined behavior since you are accessing uninitialized memory. You can easily verify this by running with the memchek module of valgrind.

Please explain why you are so set on looking for thread safety. As a general rule, LAMMPS is not thread safe. Some functionality will break if running LAMMPS in a concurrent thread. Issuing a LAMMPS command will generate an error, extract functions can be used, but may be subject to race conditions.

1 Like
3

Hi @akohlmey,

I compiled my code with the current “develop” branch from github.

I posted the above because there is no standard lammps_extract function to retrieve the interatomic potential coefficients from a LAMMPS object and I was just wondering if the line below gives me the correct sigma:

double **sigma_temp = ( double **)lmp->force->pair->extract( const_cast<char*>( "sigma" ) , 2 );

Thread-safety was probably not the correct term to use here and I apologize for that. My intention was to ask if retrieving the sigma values from lmp->force->pair->extract can lead to undefined behaviour of my code (e.g., because garbage values might be stored inside the sigma_temp array the moment I read the relevant data). However, I believe everything is clear now after reading your response, so thank you!

4

We have not added a generic function for that specific purpose because this is highly specific for different pair styles (or bond/angle/dihedral/improper). The typical way to implement access for a specific pair style is to include the header for the corresponding pair style (e.g. #include "pair_lj_cut.h")

and then have:

PairLJCut *pstyle = dynamic_cast<PairLJCut *>(force->pair);
if (!pstyle) // error out since either pair style not defined or incompatible style 
  error->all(FLERR, "Incorrect pair style");

…and from there on you can use the pstyle pointer and access its members. If needed, you make them public.

The Pair::extract() method is mostly meant for fix adapt or kspace.

Please also note, that your use of the extract function is incorrect. The second argument has to be an int variable. Its value will be set to either 0 or 1 or 2. For cutoff this matters: with 0 there is only a global cutoff with 2 there is a per pair of types cutoff.

1 Like
5

Got it! Thanks for your replies!

The suggested solution with the use of the pstyle pointer works just fine for what I am trying to do :slight_smile:!