Luke,# Okay, I figured it out. The betas in the paper you sent are the coefficients of the quintic spline function (beta0 + beta1x + beta2x^2 … beta5*x^5, six coefficients) that represents the g(cos(theta)) function. Both the Brenner 2002 and Stuart (2000) report the actual values of the g(cos(theta)) and its first and second derivatives at the spline nodes. The original codes written by Stuart and Brenner used a routine from Numerical Recipes to calculate the spline coefficients in the MD code before the simulation starts running so they did not report the coefficients of the spline function. It appears that the developer of the lammps version of the code calculated the coefficients externally and used them as input parameters. The Lindsay (Lindsay, PHYSICAL REVIEW B 81, 205441 2010) paper you sent does the same thing.
If you look at the CH.airebo file you have a section labeled “#gC1 and gC2”. Below you have the number of spline nodes, the values of cos(theta) at the nodes, and then two sets of coefficients for gc1 and gc2 (see Brenner 2002, pg 791). There are 4 sets of coefficients, with 6 coefficients in each set. In the CH.airebo file the numbers that most closely match the beta0-beta5 numbers from Lindsay are:
0.6900250000
5.4601600000
23.0108000000
54.9086400000
68.6124000000
34.7051520000
This represents the spline function for bonds between 109.5 and 120 degrees, which makes sense because Lindsay was trying to optimize graphene.
You can see that these are close to that reported in the Lindsay paper. I double checked the lammps version against my original version of the code and I get the same values. It would appear that the coefficients reported in Lindsay paper are in error.
It should be fine if you replace these numbers with the optimized numbers in the Lindsay paper.
I am not as sure about the torsional part of the potential. Again there is a spline function T that describes the interactions. This time it is a tricubic spline. It should have 4x4x4 = 64 coefficients. between each set of nodes. This is how the CH.airebo tabulates the data, so going from a single parameter, T0 to the 64 coefficients requires you to plug it into a tricubic spline. I have not done this so I can’t say how it works exactly. Don Brenner’s code: http://www.mse.ncsu.edu/CompMatSci/ will have the subroutine to calculate the tricubic spline coefficients. I think the subroutine name is tricoef or something similar. You can put some write statements in to dump out the coefficients to see how it works.
Good luck.
Dave