I want to tabulate the Finnis and Sinclair potential in LAMMPS setfl format.

As far as my knowledge, for LAMMPS, we need to provide a look-up table and LAMMPS will pick up the appropriate values for V®, rho® and F(rho) from this table by interpolation.

Now, Rho (i) = Sum over j [ ( r_ij -d ) ] ^ 2.

Rho is maximum at r = 0 and for d=4.400224 (taken from the F&S paper), the maximum value of Rho is 19.3619730879893 (Rho_Max) . So, in LAMMPS EAM table, we will tabulate F from Rho=0 to Rho= i*dRho where dRho = (Rho_max - 0)/dN (dN= number of grid points)

But, this spatial variation of Rho around site i is due to 1 atom only whose distance varies from 0 to d.

In real situation, the site i will be surrounded by many atoms (e.g. in equilibrium case, 8+6 = 14 atoms) and the density at site i will be much higher than Rho_max= 19.361973.

Now, the corresponding F(rho) is tabulated only for Rho values ranging from 0 to Rho_Max at ‘Nrho’ number of grid points where, dRho = ( Rho_Max - 0 )/Nrho.

In that case, how do LAMMPS pick up the correct value from the look up table when the effective Rho at an atomic position is more than Rho_Max ?

This isn’t a LAMMPS Q, it’s a Q about the format of EAM files

which are DYNAMO format (pre-LAMMPS). I don’t understand

your fomula with Rho as a function of r and d. It’s just a function

of r, E.g. on the LAMMPS web page: F (sum rho(rij)).

So the tabulated rho function is over r values and returns some

boudned value of rho.

The F function (embedding func) is over rho values, which are

summed over many neighbor atoms. The max value for that

sum can be anything you choose. It doens’t have to be the

rho value for one atom pair.

If you don’t think that is correct, then put build a small system

of an atom with a few neighbors, add some print statements

to pair_eam.cpp and verify that nothing is overflowing max values.