confusion about special bonds and pair potential

Hi all,

I am very confused about how “special_bonds” works and how it affects the pair potential. Say I have the following input

special_bonds fene

bond_style fene

bond_coeff 1 30.0 1.5 1.0 1.0

pair_style lj/cut 1.12

pair_modify shift yes

pair_coeff 1 1 1.0 1.0 1.12

My questions are

1,

Between bonded atoms 1 and 2, since I have bond_coeff 1 30.0 1.5 1.0 1.0 and

"pair_coeff 1 1 1.0 1.0 1.12", am I counting nonzero LJ potential for 1 and 2 twice?

If not(i.e., if there is only one nonzero LJ), which set of the coefficients are not working? **"**bond_coeff 1 30.0 1.5 1.0 1.0" or “pair_coeff 1 1 1.0 1.0 1.12”?

2,

If in the above input I do a little change

bond_coeff 1 30.0 1.5 1.0 1.0 ==> bond_coeff 1 30.0 1.5 0.0 1.0

Now between bonded atoms 1 and 2, are there still nonzero LJ potentials since I have

"special_bonds fene", **"**bond_coeff 1 30.0 1.5 0.0 1.0"

and “pair_coeff 1 1 1.0 1.0 1.12” ?

How about the nonbonded atoms, are there still nonzero LJ potential?

Thanks!!!

Hi all,

I am very confused about how "special_bonds" works and how it affects the
pair potential. Say I have the following input

special_bonds fene

bond_style fene

bond_coeff 1 30.0 1.5 1.0 1.0

pair_style lj/cut 1.12

pair_modify shift yes

pair_coeff 1 1 1.0 1.0 1.12

My questions are

1,

Between bonded atoms 1 and 2, since I have "bond_coeff 1 30.0 1.5 1.0 1.0"
and

"pair_coeff 1 1 1.0 1.0 1.12", am I counting nonzero LJ potential for 1 and
2 twice?

no.

If not(i.e., if there is only one nonzero LJ), which set of the
coefficients are not working? "bond_coeff 1 30.0 1.5 1.0 1.0" or
"pair_coeff 1 1 1.0 1.0 1.12"?

"working" is a bad choice of word. for atoms that have a (topological)
bond between them, you have an 1-2 exclusion meaning that the bonded
interaction replaced the non-bonded (i.e. pair) interaction is
excluded.

2,

If in the above input I do a little change

bond_coeff 1 30.0 1.5 1.0 1.0 ==> bond_coeff 1 30.0 1.5 0.0 1.0

Now between bonded atoms 1 and 2, are there still nonzero LJ potentials
since I have

"special_bonds fene", "bond_coeff 1 30.0 1.5 0.0 1.0"

and "pair_coeff 1 1 1.0 1.0 1.12" ?

no. you have disabled the repulsive part of the FENE bond. the
computation set up by the pair style is unaffected by this.

How about the nonbonded atoms, are there still nonzero LJ potential?

pair_style and bond_style don't know of each other's parameters.
exclusions are processed at the level of building the neighbor lists,
which encode whether a pair of atoms is a 1-2, 1-3 or 1-4 pair, i.e.
two directly bonded atoms, two atoms bonded to the same central atom,
and two atoms bonded to two atoms forming a bond, respectively.

for more details, i suggest reading through a text book that explains
how conventional classical force fields are constructed and
interactions computed.

axel.

Some more clarification

Some more clarification

1,
"you have an 1-2 exclusion meaning that the bonded interaction replaced the non-bonded (i.e. pair) interaction is excluded."
I have difficulty understanding what you really mean here. Are saying for 1-2 atoms, the interaction set by pair_style is excluded? But the lammps documentation says "In LAMMPS, pair potentials are defined between pairs of atoms that are within a cutoff distance and...".

please follow my advice and read up on how exclusion work (and why
they are used) from a proper MD text book. there is no point in me
trying to recap what is properly explained in a reference and is
something that a practitioner in MD should know and understand.

2,
"no. you have disabled the repulsive part of the FENE bond. the computation set up by the pair style is unaffected by this."

OK, the repulsive part in 1-2 atoms is disabled by "bond_coeff 1 30.0 1.5 0.0 1.0". What should I do if I want to use the pair potential only introduced by pair_style, but not the one included in the bond. E.g, for "fene/expand", an additional parameter \Delta enters the repulse part, which is not what I want. I want only the term U_fene = -0.5kR_0^2 ln(1- (r_0-\Delta)^2/R_0^2) and the normal LJ potential(without \Delta) applying to all atoms(including 1-2 atoms).
Any suggestions?

then don't use special_bonds fene, but disable all exclusions with
special_bonds lj/coul 1.0 1.0 1.0
this will leave the pairwise lj potential unchanged for *all* pairs
and then adds the attractive part of the FENE potential to all pairs
defined as a bond, provided you set the epsilon for the repulsive pair
of FENE to zero as indicated. you still have to make sure that you
choose the correct sigma, as that determines where the attractive part
starts. please note that this kind of complex manipulation may require
careful checking that you don't introduce a discontinuity in the total
potential. typically bond styles have a potential energy of zero at
the equilibrium distance, which is where FENE should switch from
repulsive to attractive.

axel.