Help with heaviside harmonic force

Dear all,

My harmonic force function for each bond is
f=k0*(r-r0) if r=rc

But I can only do

bond_style harmonic
bond_coeff * k0 r0

I wonder how to find out if any one of bond lengths r is larger than rc and change only this bond coefficients from k0 and r0 to k1 and r1, but keep others unchanged.

I would really appreciate it.

Best regards,

Alex J.

Dear all,

My harmonic force function for each bond is
f=k0*(r-r0) if r<rc
f=k1*(r-r1) if r>=rc

But I can only do
bond_style harmonic
bond_coeff * k0 r0

I wonder how to find out if any one of bond lengths r is larger than rc and
change only this bond coefficients from k0 and r0 to k1 and r1, but keep
others unchanged.

changing the functional form like you describe requires writing a new
bond style. however, you can also just use the bond_style table and
just tabulate your custom potential. that is usually the preferred
way, unless you are comfortable modifying and compiling source code. i
strongly recommend, however, to first practice using tabulation with a
bond style that *is* supported, e.g. the regular harmonic bond style
and make sure that you figured out how to generate tables, that match
the analytic function
and optimized the number of tabulation points needed and the range to
be covered.

axel.

Dear Axel,

Thank you very much for the reply.

I would like to try what you said. But would you mind telling me a little bit more about how to generate tables? I do not see the “bond_style table” anywhere and have no idea about the generation.

Best regards,

Here you go:

http://lammps.sandia.gov/doc/bond_table.html