Shake cluster are connected

Hi,

I am trying to keep my water molecules rigid in my simulation using fix. Here is the fix command

fix 1 H2O shake 0.0001 20 100 b 1 a 1

I have 1 bond type and 1 angle type. When i try to run, i am getting a error “Shake cluster are connected”

ERROR: Shake clusters are connected (…/fix_shake.cpp:1002)
Last command: fix 1 wat shake 0.0001 20 100 b 1 a 1

I am inclined to guess that you have accidentally submitted a DATA file containing additional bonds (or angles) between different water molecules. Or perhaps you have other molecules in your system which also use bond_type 1 or angle_type 1.

You should visualize your system in VMD to check that the bonds look correct.

If you cannot detect a problem, visually, then it’s time to count the number of connected sets of atoms (a.k.a. “fragments”) in your original system. THEN check to see that it matches the number of molecules that you think you have in your system. If they don’t match, then this explains why you are receiving this error message.

One way to do this is to use the “compute atom/fragment” command in LAMMPS. This command will assign each atom to a “fragment”, which you can define as any collection of atoms which are connected by bonds (of any type):

https://lammps.sandia.gov/doc/compute_cluster_atom.html

Once you have done that, you can write the list of atoms and their fragments to a file. There are several ways to do that. Then you can write your own script to count the number of disconnected fragments. Here is one way:

You could use a “set” command similar to the command below to assign unique molecule-IDs to their “fragment-ID” (calculated using compute atom/fragment), and then use a “write_data” command to create a new data file. The 2nd-column of the “Atoms” section of that data file (which normally contains the molecule-IDs), will contain the “fragment-IDs” for each atom (the fragment-number that each atom belongs to). Then count the number of unique fragment-IDs. Here is a brief sketch how to do that:

I am inclined to guess that you have accidentally submitted a DATA file
containing additional bonds (or angles) between different water molecules.
Or perhaps you have other molecules in your system which also use bond_type
1 or angle_type 1.

You should visualize your system in VMD
<https://github.com/jewettaij/moltemplate/blob/master/examples/all_atom/force_field_OPLSAA/waterSPCE%2Bmethane/README_visualize.txt>
to check that the bonds look correct.

If you cannot detect a problem, visually, then it's time to count the
number of connected sets of atoms (a.k.a. "fragments") in your original
system. THEN check to see that it matches the number of molecules that you
think you have in your system. If they don't match, then this explains why
you are receiving this error message.

One way to do this is to use the "compute atom/fragment" command in
LAMMPS. This command will assign each atom to a "fragment", which you can
define as any collection of atoms which are connected by bonds (of any
type):

https://lammps.sandia.gov/doc/compute_cluster_atom.html

Once you have done that, you can write the list of atoms and their
fragments to a file. There are several ways to do that. Then you can
write your own script to count the number of disconnected fragments. Here
is one way:

You could use a "set" command similar to the command below to assign
unique molecule-IDs to their "fragment-ID" (calculated using compute
atom/fragment), and then use a "write_data" command to create a new data
file. The 2nd-column of the "Atoms" section of that data file (which
normally contains the molecule-IDs), will contain the "fragment-IDs" for
each atom (the fragment-number that each atom belongs to). Then count the
number of unique fragment-IDs. Here is a brief sketch how to do that:

---------------

read_data "orig_file.data"
compute CWhichMol all fragment/atom
set mol c_CWhichMol
run 0 # <-- this might be optional
write_data "new_mIDs.data"

---------------

Warning: I have not tested these commands. For details, see the
documentation at:
https://lammps.sandia.gov/doc/set.html
https://lammps.sandia.gov/doc/write_data.html

Then extract the "Atoms" section of the data file you just created (eg
"new_mIDs.data"), and count the number of unique integers which appear in
the 2nd column of that section of the file (which is where the molecule-ID
numbers are stored, assuming your are using "atom_style full").

---------------

extract_lammps_data.py Atoms < new_mIDs.data > new_Atoms_section.txt
sort -g -k 3 < new_Atoms_section.txt
    > awk 'BEGIN{i=-1; n=0} {if ($2>i) {n++} i=$2} END{print n}'

I'm sorry. Typo. That should have been:

sort -g -k 2 < new_Atoms_section.txt
    > awk 'BEGIN{i=-1; n=0} {if ($2>i) {n++} i=$2} END{print n}'

There are probably other mistakes in these instructions as well. Be
careful.

Hi Andrew,

Thanks for your response. You were absolutely correct. My data file had some H2O having more coordination. Your scripts helped me to figure out it. Thanks alot for it.

Cheers,
Naresh