Insufficient Jacobi rotations for rigid body

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

I just run a very simple simulation, with 27 water molecules. It gives the error: Insufficient Jacobi rotations for rigid body. I have no idea how to solve this problem. Thanks for your time and help!

Regards,
Joy

Here is the In file:

units real
atom_style full
dimension 3
boundary p p p

read_data data.water

velocity all create 300 432567 dist uniform

group water molecule <> 1 27

pair_style lj/cut/coul/long 10.0
kspace_style ewald 1.0e-4

pair_coeff 1 1 0.2106 3.166
pair_coeff * 2 0.0000 0.0000
pair_coeff * 3 0.0000 0.0000

fix 1 water rigid/nvt molecule temp 300.0 300.0 100.0
fix 2 all temp/rescale 50 300 300 50 1.0

neighbor 2.0 bin
neigh_modify exclude type 1 2
neigh_modify exclude type 1 3
neigh_modify delay 0 every 1 one 2000 page 20000

timestep 0.1

thermo_style custom step temp epair ebond eangle vol
thermo 500
dump myDump all xyz 1000 dump.dat
run 50000

write_restart water.restart

And data file:

LAMMPS Atom File

108 atoms
0 bonds
0 angles
0 dihedrals
0 impropers

3 atom types
0 bond types
0 angle types

-4.72 4.72 xlo xhi
-4.72 4.72 ylo yhi
-4.72 4.72 zlo zhi

Masses

1 15.9994
2 1.008
3 1e-10

Atoms

1 1 1 0.000 -3.147 -3.147 -3.147
2 1 2 0.410 -2.458 -3.057 -3.806
3 1 2 0.410 -2.684 -3.112 -2.310
4 1 3 -0.820 -2.992 -3.130 -3.123
5 2 1 0.000 -0.000 -3.147 -3.147
6 2 2 0.410 -0.183 -2.581 -3.897
7 2 2 0.410 0.928 -3.006 -2.961
8 2 3 -0.820 0.100 -3.052 -3.223
9 3 1 0.000 3.147 -3.147 -3.147
10 3 2 0.410 3.810 -2.980 -3.816
11 3 2 0.410 3.647 -3.398 -2.370
12 3 3 -0.820 3.303 -3.158 -3.132
13 4 1 0.000 -3.147 -0.000 -3.147
14 4 2 0.410 -2.717 0.604 -3.752
15 4 2 0.410 -2.748 0.188 -2.297
16 4 3 -0.820 -3.035 0.107 -3.114
17 5 1 0.000 -0.000 -0.000 -3.147
18 5 2 0.410 0.536 -0.165 -3.923
19 5 2 0.410 0.570 -0.222 -2.410
20 5 3 -0.820 0.149 -0.052 -3.152
21 6 1 0.000 3.147 -0.000 -3.147
22 6 2 0.410 3.332 0.508 -3.936
23 6 2 0.410 3.310 -0.908 -3.402
24 6 3 -0.820 3.194 -0.054 -3.287
25 7 1 0.000 -3.147 3.147 -3.147
26 7 2 0.410 -2.554 3.382 -3.860
27 7 2 0.410 -2.724 3.492 -2.360
28 7 3 -0.820 -3.010 3.225 -3.137
29 8 1 0.000 -0.000 3.147 -3.147
30 8 2 0.410 0.413 3.812 -3.698
31 8 2 0.410 0.412 3.256 -2.290
32 8 3 -0.820 0.111 3.251 -3.106
33 9 1 0.000 3.147 3.147 -3.147
34 9 2 0.410 4.057 3.217 -3.434
35 9 2 0.410 3.155 3.463 -2.243
36 9 3 -0.820 3.270 3.199 -3.064
37 10 1 0.000 -3.147 -3.147 -0.000
38 10 2 0.410 -3.078 -2.614 -0.792
39 10 2 0.410 -2.371 -2.917 0.512
40 10 3 -0.820 -3.033 -3.044 -0.038
41 11 1 0.000 -0.000 -3.147 -0.000
42 11 2 0.410 0.509 -2.535 -0.532
43 11 2 0.410 0.202 -2.904 0.904
44 11 3 -0.820 0.096 -3.032 0.050
45 12 1 0.000 3.147 -3.147 -0.000
46 12 2 0.410 3.279 -3.086 -0.946
47 12 2 0.410 3.722 -3.860 0.277
48 12 3 -0.820 3.242 -3.234 -0.090
49 13 1 0.000 -3.147 -0.000 -0.000
50 13 2 0.410 -2.542 0.132 -0.730
51 13 2 0.410 -2.582 -0.068 0.770
52 13 3 -0.820 -2.989 0.009 0.005
53 14 1 0.000 -0.000 -0.000 -0.000
54 14 2 0.410 0.467 -0.420 -0.723
55 14 2 0.410 0.125 -0.592 0.742
56 14 3 -0.820 0.080 -0.136 0.003
57 15 1 0.000 3.147 -0.000 -0.000
58 15 2 0.410 3.669 0.489 -0.636
59 15 2 0.410 3.627 -0.819 0.126
60 15 3 -0.820 3.282 -0.044 -0.069
61 16 1 0.000 -3.147 3.147 -0.000
62 16 2 0.410 -3.124 3.312 -0.943
63 16 2 0.410 -2.298 2.750 0.195
64 16 3 -0.820 -3.029 3.116 -0.101
65 17 1 0.000 -0.000 3.147 -0.000
66 17 2 0.410 0.009 3.365 -0.932
67 17 2 0.410 0.808 3.528 0.344
68 17 3 -0.820 0.110 3.227 -0.079
69 18 1 0.000 3.147 3.147 -0.000
70 18 2 0.410 3.627 3.290 -0.815
71 18 2 0.410 3.827 3.044 0.665
72 18 3 -0.820 3.303 3.152 -0.020
73 19 1 0.000 -3.147 -3.147 3.147
74 19 2 0.410 -2.868 -2.324 2.745
75 19 2 0.410 -2.811 -3.102 4.042
76 19 3 -0.820 -3.064 -3.030 3.213
77 20 1 0.000 -0.000 -3.147 3.147
78 20 2 0.410 0.137 -2.505 2.450
79 20 2 0.410 0.719 -2.993 3.759
80 20 3 -0.820 0.115 -3.040 3.135
81 21 1 0.000 3.147 -3.147 3.147
82 21 2 0.410 3.292 -3.435 2.246
83 21 2 0.410 4.014 -2.887 3.458
84 21 3 -0.820 3.283 -3.151 3.067
85 22 1 0.000 -3.147 -0.000 3.147
86 22 2 0.410 -2.401 0.357 2.664
87 22 2 0.410 -3.199 0.536 3.938
88 22 3 -0.820 -3.053 0.120 3.188
89 23 1 0.000 -0.000 -0.000 3.147
90 23 2 0.410 0.359 -0.592 2.486
91 23 2 0.410 0.298 -0.362 3.981
92 23 3 -0.820 0.088 -0.128 3.170
93 24 1 0.000 3.147 -0.000 3.147
94 24 2 0.410 3.699 0.100 2.372
95 24 2 0.410 3.760 -0.194 3.855
96 24 3 -0.820 3.304 -0.013 3.138
97 25 1 0.000 -3.147 3.147 3.147
98 25 2 0.410 -2.637 2.917 2.370
99 25 2 0.410 -2.507 3.160 3.858
100 25 3 -0.820 -2.992 3.118 3.138
101 26 1 0.000 -0.000 3.147 3.147
102 26 2 0.410 0.807 3.211 2.637
103 26 2 0.410 0.294 2.979 4.042
104 26 3 -0.820 0.148 3.133 3.199
105 27 1 0.000 3.147 3.147 3.147
106 27 2 0.410 3.964 3.249 2.659
107 27 2 0.410 3.418 2.809 4.000
108 27 3 -0.820 3.293 3.115 3.196

I got the output as below:

Scanning data file …
Reading data file …
orthogonal box = (-4.72 -4.72 -4.72) to (4.72 4.72 4.72)
1 by 1 by 1 processor grid
108 atoms
Finding 1-2 1-3 1-4 neighbors …
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
108 atoms in group water
27 rigid bodies with 108 atoms
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62
Setting up run …
Memory usage per processor = 2.09557 Mbytes
Step Temp E_pair E_bond E_angle Volume
0 178.16984 -5961.4573 0 0 841.23238
1000 300 -6051.303 0 0 841.23238
2000 300 -6063.6642 0 0 841.23238
3000 300 -6051.1059 0 0 841.23238
4000 300 -6059.3937 0 0 841.23238
5000 300 -6037.6275 0 0 841.23238
6000 300 -6003.0647 0 0 841.23238
7000 300 -6013.4112 0 0 841.23238
8000 300 -5925.6174 0 0 841.23238
9000 300 -5863.6547 0 0 841.23238
10000 300 -5534.0281 0 0 841.23238
11000 300 -5695.5609 0 0 841.23238
12000 300 -4890.5809 0 0 841.23238
13000 300 -4897.2194 0 0 841.23238
14000 300 -4274.6457 0 0 841.23238
14500 300 280.51906 0 0 841.23238
15000 nan nan 0 0 841.23238
15500 nan nan 0 0 841.23238
16000 nan nan 0 0 841.23238
…

Loop time of 111.471 on 1 procs for 50000 steps with 108 atoms

Pair time () = 88.1273 (79.0585) Bond time () = 0.00457811 (0.004107)
Kspce time () = 4.22416 (3.78947) Neigh time () = 14.3117 (12.839)
Comm time () = 2.81155 (2.52223) Outpt time () = 0.00582933 (0.00522946)
Other time (%) = 1.98588 (1.78152)

Nlocal: 108 ave 108 max 108 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4550 ave 4550 max 4550 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 31188 ave 31188 max 31188 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 31188
Ave neighs/atom = 288.778
Ave special neighs/atom = 0
Neighbor list builds = 5684
Dangerous builds = 3478
System init for write_restart …
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62
ERROR: Insufficient Jacobi rotations for rigid body

2

Hi all,

joy,

it is very nice that you provide a complete input
to try reproducing your problem, but when reporting
something like this, please also always let us
know which lammps version you are using. thanks.

I just run a very simple simulation, with 27 water molecules. It gives the
error: Insufficient Jacobi rotations for rigid body. I have no idea how to

the problem lies elsewhere. didn't you notice all the NaN outputs.
those are an indication, that something "really bad(tm)"
has happened during the calculation of the forces.

solve this problem. Thanks for your time and help!

Regards,
Joy

Here is the In file:

units real
atom_style full
dimension 3
boundary p p p

read_data data.water

velocity all create 300 432567 dist uniform

group water molecule <> 1 27

pair_style lj/cut/coul/long 10.0
kspace_style ewald 1.0e-4

pair_coeff 1 1 0.2106 3.166
pair_coeff * 2 0.0000 0.0000
pair_coeff * 3 0.0000 0.0000

fix 1 water rigid/nvt molecule temp 300.0 300.0 100.0

ok. so far.

fix 2 all temp/rescale 50 300 300 50 1.0

ouch! this is a problem. you should not modify velocities twice.
fix rigid/nvt already contains a thermostat algorithm.
on top of that, temp/rescale is a horiible way to "control"
temperature and should only be used under extreme
circumstances, when there is no alternative. just
remove that line from your input.

neighbor 2.0 bin
neigh_modify exclude type 1 2
neigh_modify exclude type 1 3

you may also want to exclude the intramolecular
forces via the following:

neigh_modify exclude molecule water

neigh_modify delay 0 every 1 one 2000 page 20000

why do you mess with the neighbor list table settings?
the defaults for one and page should work fine. since
rigid particles move somewhat slower, you can tweak
the settings to check for rebuilds less often. e.g.
try: neigh_modify delay 4 every 2

timestep 0.1

this is a _very_ safe choice for a time step.
you can probably crank this up to 0.5fs
without much of a problem and 5x as much
work done in the same time. :wink:

thermo_style custom step temp epair ebond eangle vol
thermo 500
dump myDump all xyz 1000 dump.dat
run 50000

you input runs fine for me with and without the
modifications that i am suggesting.

what platform/compiler and version of LAMMPS are you using?

perhaps you get hit by a bug that has already been fixed.

axel.

3

Hello Axel,

Thanks very much for the detailed explanation. I revised my input as you suggested. Without “fix 2 all temp/rescale 50 300 300 50 1.0”, I still got the same error. And also, the temperature went out of control. See below:

Step Temp E_pair E_bond E_angle Volume
0 178.16984 -1409.0583 0 0 841.23238
500 603162.8 -1477.9058 0 0 841.23238
1000 596899.65 -1495.4018 0 0 841.23238
1500 550785.75 -1505.0602 0 0 841.23238
5000 391481.09 -1491.8248 0 0 841.23238
8000 680259.67 -1388.7492 0 0 841.23238
10000 968883.14 -1305.8431 0 0 841.23238
11500 51920304 -882.05436 0 0 841.23238
12000 1.918702e+08 -536.76259 0 0 841.23238
14000 1.2174756e+09 1235.4643 0 0 841.23238
15000 3.7109656e+09 2875.333 0 0 841.23238
16000 nan nan 0 0 841.23238
16500 nan nan 0 0 841.23238
17000 nan nan 0 0 841.23238
17500 nan nan 0 0 841.23238
18000 nan nan 0 0 841.23238
…

Total # of neighbors = 30157
Ave neighs/atom = 279.231
Ave special neighs/atom = 0
Neighbor list builds = 3135
Dangerous builds = 1808
System init for write_restart …
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62
ERROR: Insufficient Jacobi rotations for rigid body

By the way, I use the version of LAMMPS (4 Jul 2010). It is installed in the cluster of our department. I can not compile the program myself.

LOL.

Then I tried the same input and data file in my own computer (Ubuntu 11.04 with LAMMPS (15 Jan 2012)), everything is fine till 50000 steps. I attached the output for your interest.

Step Temp E_pair E_bond E_angle Volume
0 178.16984 -1409.0583 0 0 841.23238
500 384.23404 -1468.0202 0 0 841.23238
1000 406.93099 -1485.7684 0 0 841.23238
1500 429.32135 -1497.7881 0 0 841.23238
2000 383.44793 -1502.109 0 0 841.23238

48000 290.93562 -1503.4614 0 0 841.23238
48500 313.92916 -1503.7504 0 0 841.23238
49000 332.8507 -1502.9425 0 0 841.23238
49500 338.92999 -1501.0475 0 0 841.23238
50000 398.77662 -1508.1372 0 0 841.23238
Loop time of 148.629 on 2 procs for 50000 steps with 108 atoms

Pair time () = 30.0059 (20.1884) Bond time () = 0.0117188 (0.00788457)
Kspce time () = 28.9883 (19.5038) Neigh time () = 0.339844 (0.228653)
Comm time () = 80.457 (54.1328) Outpt time () = 0.0546875 (0.0367947)
Other time (%) = 8.77148 (5.9016)

Nlocal: 54 ave 57 max 51 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Nghost: 3910.5 ave 3918 max 3903 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Neighs: 15183.5 ave 16958 max 13409 min
Histogram: 1 0 0 0 0 0 0 0 0 1

Total # of neighbors = 30367
Ave neighs/atom = 281.176
Ave special neighs/atom = 0
Neighbor list builds = 178
Dangerous builds = 0
System init for write_restart …
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62

I still need to work on the cluster when running very big systems :(.

Another quick question:
INPUT 1:
group water molecule <> 1 27
fix 1 water rigid/nvt molecule temp 300.0 300.0 100.0

INPUT 2:

group water1 id <> 1 4
group water2 id <> 5 8
…
group water26 id <> 101 104
group water27 id <> 105 108

fix 1 water1 rigid/nvt molecule temp 300.0 300.0 100.0
fix 2 water2 rigid/nvt molecule temp 300.0 300.0 100.0
…
fix 26 water26 rigid/nvt molecule temp 300.0 300.0 100.0
fix 27 water27 rigid/nvt molecule temp 300.0 300.0 100.0

Are they the same setup? I learn INPUT 2 from the example given in the rigid folder in LAMMPS.

Thanks again for your help.

Joy

4

Hello Axel,

Thanks very much for the detailed explanation. I revised my input as you suggested. Without “fix 2 all temp/rescale 50 300 300 50 1.0”, I still got the same error. And also, the temperature went out of control. See below:

Step Temp E_pair E_bond E_angle Volume
0 178.16984 -1409.0583 0 0 841.23238
500 603162.8 -1477.9058 0 0 841.23238
1000 596899.65 -1495.4018 0 0 841.23238
1500 550785.75 -1505.0602 0 0 841.23238
5000 391481.09 -1491.8248 0 0 841.23238
8000 680259.67 -1388.7492 0 0 841.23238
10000 968883.14 -1305.8431 0 0 841.23238
11500 51920304 -882.05436 0 0 841.23238
12000 1.918702e+08 -536.76259 0 0 841.23238
14000 1.2174756e+09 1235.4643 0 0 841.23238
15000 3.7109656e+09 2875.333 0 0 841.23238
16000 nan nan 0 0 841.23238
16500 nan nan 0 0 841.23238
17000 nan nan 0 0 841.23238
17500 nan nan 0 0 841.23238
18000 nan nan 0 0 841.23238
…

Total # of neighbors = 30157
Ave neighs/atom = 279.231
Ave special neighs/atom = 0
Neighbor list builds = 3135
Dangerous builds = 1808
System init for write_restart …
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62
ERROR: Insufficient Jacobi rotations for rigid body

By the way, I use the version of LAMMPS (4 Jul 2010). It is installed in the cluster of our department. I can not compile the program myself.

Well, then you either better learn how to compile lammps yourself or make somebody do it for you. That old version is obviously having a bug.

LOL.

Then I tried the same input and data file in my own computer (Ubuntu 11.04 with LAMMPS (15 Jan 2012)), everything is fine till 50000 steps. I attached the output for your interest.

Step Temp E_pair E_bond E_angle Volume
0 178.16984 -1409.0583 0 0 841.23238
500 384.23404 -1468.0202 0 0 841.23238
1000 406.93099 -1485.7684 0 0 841.23238
1500 429.32135 -1497.7881 0 0 841.23238
2000 383.44793 -1502.109 0 0 841.23238

48000 290.93562 -1503.4614 0 0 841.23238
48500 313.92916 -1503.7504 0 0 841.23238
49000 332.8507 -1502.9425 0 0 841.23238
49500 338.92999 -1501.0475 0 0 841.23238
50000 398.77662 -1508.1372 0 0 841.23238
Loop time of 148.629 on 2 procs for 50000 steps with 108 atoms

Pair time () = 30.0059 (20.1884) Bond time () = 0.0117188 (0.00788457)
Kspce time () = 28.9883 (19.5038) Neigh time () = 0.339844 (0.228653)
Comm time () = 80.457 (54.1328) Outpt time () = 0.0546875 (0.0367947)
Other time (%) = 8.77148 (5.9016)

Nlocal: 54 ave 57 max 51 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Nghost: 3910.5 ave 3918 max 3903 min
Histogram: 1 0 0 0 0 0 0 0 0 1
Neighs: 15183.5 ave 16958 max 13409 min
Histogram: 1 0 0 0 0 0 0 0 0 1

Total # of neighbors = 30367
Ave neighs/atom = 281.176
Ave special neighs/atom = 0
Neighbor list builds = 178
Dangerous builds = 0
System init for write_restart …
Ewald initialization …
G vector = 0.273155
vectors: actual 1d max = 40 2 62

I still need to work on the cluster when running very big systems :(.

Another quick question:
INPUT 1:
group water molecule <> 1 27
fix 1 water rigid/nvt molecule temp 300.0 300.0 100.0

INPUT 2:

group water1 id <> 1 4
group water2 id <> 5 8
…
group water26 id <> 101 104
group water27 id <> 105 108

fix 1 water1 rigid/nvt molecule temp 300.0 300.0 100.0
fix 2 water2 rigid/nvt molecule temp 300.0 300.0 100.0
…
fix 26 water26 rigid/nvt molecule temp 300.0 300.0 100.0
fix 27 water27 rigid/nvt molecule temp 300.0 300.0 100.0

Are they the same setup? I learn INPUT 2 from the example given in the rigid folder in LAMMPS.

No, they are not the same. The water molecules are coupled to thermostats differently, so your dynamics will be different.

Axel