[lammps-users] Using FENE potential to constraint motion


I wish to setup a wall with atoms vibrating about their mean positions constrained by a fene potential. For the same I am creating a molecule of 2 atoms, connected by a fene bond and keeping one atom frozen at the expected mean position. To avoid any other pair potential interaction I am setting pair_coeff 1 2 0.0 0.0. But upon simulating, the program gives warnings that FENE bond too long. I am unable to find out the reason for this, as the atoms are interacting only through the FENE potential and there should be no reason that they are repelled out so far. My input file is as follows. I would be extremely thankful if somebody can point me to the idea I might be missing out here.

3-d LJ flow simulation

units lj
dimension 3
newton on
processors 2 2 1
boundary p p p
atom_style bond

create geometry

read_data data.Capillary
special_bonds fene
bond_style fene
bond_coeff 1 30 1.5 0.0 0.0

define groups

group fix type 1 # frozen atom
group mobile type 2 #mobile atom
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
neigh_modify exclude group fix fix check no #exclude all interactions between frozen atoms

LJ potentials

pair_style lj/cut 2.5
pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.8 0.6

initial velocities

communicate single vel yes
velocity mobile create 1.0 1290
velocity fix create 0.0 1290
fix 1 mobile nvt 1.0 1.0 100.0
fix 2 fix setforce 0.0 0.0 0.0


timestep 0.00144
run_style verlet
thermo 100
reset_timestep 0
dump 1 all xyz 500 Capillary.xyz
run 10000

With regards,

Manik Mayur
Graduate student
Microfluidics Lab
Dept. of Mechanical Engg.
IIT Kharagpur

bond_style fene
bond_coeff 1 30 1.5 0.0 0.0

What you describe should be possible. I've done
it myself. I would start with a simple system, monitor
the bond lenghts, viz what is happening, etc.

The FENE coeff command above is not typical. You
need epsilon,sigma settings to get the bond length
right. You could also try to tether the wall atoms
to the frozen ones with bond_style harmonic, with
a r0 = 0.0.


Hi Steve,

Thanks a lot for your reply. I used 0.0 as epsilon and sigma just to make sure that the dummy atom and the actual atom do not have a LJ interaction and just the FENE interaction. If I give them some value other than zero, the system blows up with inf. potential energy and forces. Also, I was not able to understand why the FENE bonds are extending beyond the limit if it is enforced not to. I do not want to use harmonic bonds in this problem as I am trying to reproduce some published results which use FENE.

If you have used this idea earlier, I would be grateful if you send me that input file if it is ok with you.

Thanks and regards,

Manik Mayur
Graduate student
Microfluidics Lab
Dept. of Mechanical Engg.
IIT Kharagpur

2009/12/23 Steve Plimpton <[email protected]>

I don't have the input file anymore. I would try it with harmonic
bonds. But there is no reason FENE shouldn't work - you just
need to make sure the bonds don't stretch too much. Unlike
harmonic bonds, FENE bonds have a singularity at maximum