[lammps-users] The fix shake command

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I am using fix shake command for CH2 group, cyclohexane (CH2)6. But I got the error message

“Shake clusters are connected”

I checked this means “A single cluster specified by the fix shake command must have a single central atom with up to 3 other atoms bonded to it”. However, in this specific ring structure molecule, how can be only one central atom? Can anyone tell me how should I use this command for cyclohexane?

And another question, is there any command in LAMMPS can exclude the intramolecular pair-wise interactions? I only want angle_style and dihedral_style potentials are used within one molecule, and the bond length is fixed.

Your help would be most appreciated.

Best
Yajie

2

You can't do general SHAKE with LAMMPS, e.g. of
large moleclues. There is a fix rigid, which you can
look at to see if might do what you want. Look at
the neigh_modify exclude command for your 2nd Q.

Steve

3

Dear Steve,

Thank you very much for your reply.

Yes, the neigh_modify does help for the exclusion of intramolecular pair-wise interaction. But the fix rigid may not work for what I want, because I only want the bond length is fixed, but there still are bond bending and angle torsion within one molecule. So I am wondering can LAMMPS do this?

Most appreciate

Best
Yajie

4

Dear Steve,

Thank you very much for your reply.

Yes, the neigh_modify does help for the exclusion of intramolecular
pair-wise interaction. But the fix rigid may not work for what I want,
because I only want the bond length is fixed, but there still are bond
bending and angle torsion within one molecule. So I am wondering can LAMMPS
do this?

yajie,

how "fixed" do you need your bond lengths? a simple way to "emulate"
constraint is to crank up the force constant. for that you probably have
to pay a price by reducing the time step, but you can make up for it
by using the respa integrator. for example similar to this:

in one of my simulations i have replaced:

timestep 5.0
run_style verlet
neigh_modify delay 5

with:
timestep 20.0
run_style respa 2 4 bond 1 angle 1 pair 2
neigh_modify delay 1

and can run systems with pretty stiff bonds and
relatively soft non-bonded potentials very efficiently.

cheers,
    axel.

5

Thank you very much, Axel.

I am referring the paper Journal of Chemical Physics, 111, 9731-9738, 1999 and using its potential for cyclohexane. Within one molecule, the bond length is mentioned to be fixed at 1.535A. And the model consists of variable bond bending angles and torsion angles. The pair-wise potential is only considered for nonbonded interactions between different molecules. And this pair-wise potential is indeed infinite huge at small separation (about less than 3.5A). So I just follow the model and use LAMMPS to work on it.

I have excluded the intramolecular pair-wise potential by using neigh_modify and bond_style with none, but I am not sure whether it can realize what I want. I still have bond_style with harmonic and dihedral_style with opls for bond bending and torsion respectively.

Do you think my case can be also worked out with your example? If you have any more suggestion, could you please tell me? Most appreciate.

Best
Yajie

6

Thank you very much, Axel.

I am referring the paper Journal of Chemical Physics, 111, 9731-9738,
1999 and using its potential for cyclohexane. Within one molecule, the
bond length is mentioned to be fixed at 1.535A. And the model consists

i don't think it would make much of a difference, if you just use
a regular C-C harmonic force constant and see, if you get a reasonable
average distance. if not, just crank up the force constant by a
factor of 10 and see what happens then.

i don't have the time to look up the potential, but i assume it
is using a united atom approach, i.e. there is only one site
for the carbon and no explicit hydrogens, right?

of variable bond bending angles and torsion angles. The pair-wise
potential is only considered for nonbonded interactions between
different molecules. And this pair-wise potential is indeed infinite

yes. that is normal. to handle this, all you have to do is to set
special_bonds lj/coul 0.0 0.0 0.0

huge at small separation (about less than 3.5A). So I just follow the
model and use LAMMPS to work on it.

I have excluded the intramolecular pair-wise potential by using
neigh_modify and bond_style with none, but I am not sure whether it

if you have a united atom model, the special bonds command from
above should already work for you in cyclohexane, as the longest
atom-atom distance is 1-4.

can realize what I want. I still have bond_style with harmonic and
dihedral_style with opls for bond bending and torsion respectively.

if you define the bond as harmonic and use an arbitrary force constant,
say 10 kcal/mol and then use fix shake on the cyclohexane molecule
but only on this one bond type it should work, if i understand
the documentation correctly.

Do you think my case can be also worked out with your example? If you
have any more suggestion, could you please tell me? Most appreciate.

i think both variants should work reasonably well. i suspect you are
just confused about how these different options all play together.

i suggest you set up a system with just a few molecules and try out
the different variants, _systematically_. if you do that, i am certain,
you can find a way that works and that is sufficiently accurate.

cheers,
   axel.

7

Thank you very much, Axel.

Yes, it is a united atom for one CH2 group and one cyclohexane molecule has 6 CH2 groups.

I did try fix shake on the cyclohexane molecule, but it didn’t work because a single cluster specified by the fix shake command must have a single central atom with up to 3 other atoms bonded to it. So I can’t do general SHAKE with LAMMPS, e.g. of large moleclues as Steve mentioned in the earlier reply.

Anyway, I will try as much as I can. Thanks again for your suggestion.

Best
Yajie

8

Thank you very much, Axel.

Yes, it is a united atom for one CH2 group and one cyclohexane
molecule has 6 CH2 groups.

I did try fix shake on the cyclohexane molecule, but it didn't work
because a single cluster specified by the fix shake command must have
a single central atom with up to 3 other atoms bonded to it. So I
can't do general SHAKE with LAMMPS, e.g. of large moleclues as Steve
mentioned in the earlier reply.

ok. i didn't need to use SHAKE in LAMMPS before, so i didn't notice.

anyways with the special_bonds command i mentioned and a force
constant of 25-50.0 kcal/mol (or more) you should get fairly stiff
CH2-CH2 bonds and no pairwise interactions within each molecule,
but keep the angles and dihedral interactions.

i would first use a short time step and see how well it works and
ramp up the time step until you are running into energy conservation
or stability problems. go a step back, and then switch to respa
and double the outer time step and check again and so on.

cheers,
   axel.