Dear LAMMPS developers,

I’m treating a polydisperse system for which a new pair and bond styles are needed.

For this, I’ve modified the lj/cut and the fene potentials (both the .cpp and .h files) in order to add another contribution that depend on a second cutoff. The entire pair and bond potentials should be the following:

pair_style

V_WCA=4*epsilon*((sigma/r)^12-(sigma/r)^6)+epsilon if r<2^(1/6)*sigma

V_alpha=-epsilon*alpha if r<2^(1/6)*sigma

V_alpha=0.5*epsilon*alpha*(cos(gamma*(r/sigma)^2+beta)-1) if 2^(1/6)*sigma<r<R_0*sigma

bond_style

V_fene+V_alpha where V_alpha is defined above

where gamma and beta are constants, epsilon is fixed to 1.0 directly in the source code, alpha is provided in the input script instead of epsilon and V_fene is the one already built-in in lammps.

To introduce a second cutoff (R_0*sigma), I’ve multiplied cutsq in the source code by a factor FACT=1.5*1.5/2^(1/3) (I tried also the other way around, i.e. to provide R_0*sigma in the pair_coeff command and get the other cutoff as 2^(1/3)*1.5*1.5).

Apart from substituting appropriately epsilon with 1.0 and alpha with epsilon in the source code, I’ve modified the compute function of the pair and bond styles. Lammps then compiles correctly without errors.

However, when starting the simulation, "Non-numeric atom coords - simulation unstable” appears and from tests I’ve performed, it appears that the problem resides in the pair_style.

What are the reasons why this shouldn’t work? Am I missing any step in the modification?

I attach the compute functions I’ve modified.

In the input script:

bond_style fene/alpha

special_bonds fene

pair_style lj/cut/alpha 1.122462048309

pair_modify shift yes

pair_modify mix arithmetic

read_data init.dat

where bond and pair coefficients are provided in the init.dat file.

Thank you.

Best regards,

Fabrizio

computefunctions.cpp (5.82 KB)