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)