[lammps-users] Ambiguity in lammps

_Sabine_Leroch · August 28, 2009, 12:25pm

Dear lammps users

Up to now I created data files with msi2lpm using .car and .mdf files. Now I switched to the ch2lmp tool using pdb and psf data. Although the molecule geometries, force parameters, bonds, angles, box sizes, surrounding water molecules etc. are exactely the same in the two data files (only the atom, typeIDs etc. are differently distributed among the atoms), there are small differences in the resulting structures.

So I looked at the log.files and the strange thing about it is :although I used the same in. files, Lammps versions and machine its header is not the same (have a look below). What is most confusing the max numbers of bonds and dihedrals are to low in the ch2lmp log.file. The calculated numbers (for instance max bonds/atom=3) cant be true as can be seen by comparison with the Bonds list in the ch2lpm data file (for instance atomId 10 has 4 next neighbors).

What is going on here? Why differences in the logfiles? Has the difference in max bond consequences on the results, or is it just an information not used any more later?

Hope somebody can help me.

Best regards

Sabine

head of log.file for the msi datafile:

LAMMPS (9 Jan 2009)

Scanning data file …

4 = max bonds/atom

6 = max angles/atom

15 = max dihedrals/atom

2 = max impropers/atom

Reading data file …

orthogonal box = (0 -10 0) to (32 32 32.86)

from data.file:

Bonds

1 1 1 5 # CT1 NH3

2 2 1 2 # HC NH3

3 2 1 3 # HC NH3

4 2 1 4 # HC NH3

5 3 5 6 # CT1 HB

6 4 5 7 # CC CT1

7 5 5 8 # CT1 CT2

8 6 7 23 # CC OC

9 6 7 22 # CC OC

10 7 8 9 # CT2 HA

11 7 8 10 # CT2 HA

12 8 8 11 # CA CT2

13 9 11 12 # CA CA

14 9 11 20 # CA CA

15 10 12 13 # CA HP

16 9 12 14 # CA CA

17 10 14 15 # CA HP

18 9 14 16 # CA CA

19 10 16 17 # CA HP

20 9 16 18 # CA CA

21 10 18 19 # CA HP

22 9 18 20 # CA CA

23 10 20 21 # CA HP

head for the ch2lmp data.file:

LAMMPS (9 Jan 2009)

units real

atom_style full

boundary p p f

pair_style lj/charmm/coul/long 15.0 16.0 16.0

bond_style harmonic

angle_style charmm

dihedral_style charmm

improper_style harmonic ## this block is missing in the above version

kspace_style pppm 0.0001

kspace_modify slab 3.0

read_data phe_car.data # here is the name of the file above just sanning data file

3 = max bonds/atom # wrong input !!!

6 = max angles/atom

12 = max dihedrals/atom # wrong input

2 = max impropers/atom

orthogonal box = (0 -10 0) to (32 32 32.86)

data.file:

Bonds

1 5 1 3 # CC OC

2 5 1 2 # CC OC

3 4 1 8 # CC CT1

4 8 4 8 # CT1 NH3

5 10 5 4 # HC NH3

6 10 6 4 # HC NH3

7 10 7 4 # HC NH3

8 7 8 9 # CT1 HB

9 6 10 8 # CT1 CT2

10 9 10 11 # CT2 HA

11 9 10 12 # CT2 HA

12 2 13 10 # CA CT2

13 3 14 15 # CA HP

14 1 14 13 # CA CA

15 1 16 14 # CA CA

16 3 16 17 # CA HP

17 1 18 22 # CA CA

18 1 18 16 # CA CA

19 3 18 19 # CA HP

20 1 20 13 # CA CA

21 3 20 21 # CA HP

22 1 22 20 # CA CA

23 3 22 23 # CA HP

_Sabine_Leroch · August 28, 2009, 2:39pm

Hi

I found the first difference by myself. The first block in the msi log.file is not missing. Just copied the log.file entry reading from restart not from data. Sorry

But the discrepancy concerning the max bonds/atoms remains.

Greetings

Sabine

“Sabine Leroch” <sabine.leroch@…1181…> 8/28/2009 2:25 PM >>>

Dear lammps users

Up to now I created data files with msi2lpm using .car and .mdf files. Now I switched to the ch2lmp tool using pdb and psf data. Although the molecule geometries, force parameters, bonds, angles, box sizes, surrounding water molecules etc. are exactely the same in the two data files (only the atom, typeIDs etc. are differently distributed among the atoms), there are small differences in the resulting structures.

So I looked at the log.files and the strange thing about it is :although I used the same in. files, Lammps versions and machine its header is not the same (have a look below). What is most confusing the max numbers of bonds and dihedrals are to low in the ch2lmp log.file. The calculated numbers (for instance max bonds/atom=3) cant be true as can be seen by comparison with the Bonds list in the ch2lpm data file (for instance atomId 10 has 4 next neighbors).

What is going on here? Why differences in the logfiles? Has the difference in max bond consequences on the results, or is it just an information not used any more later?

Hope somebody can help me.

Best regards

Sabine

head of log.file for the msi datafile:

LAMMPS (9 Jan 2009)

Scanning data file …

4 = max bonds/atom

6 = max angles/atom

15 = max dihedrals/atom

2 = max impropers/atom

Reading data file …

orthogonal box = (0 -10 0) to (32 32 32.86)

from data.file:

Bonds

1 1 1 5 # CT1 NH3

2 2 1 2 # HC NH3

3 2 1 3 # HC NH3

4 2 1 4 # HC NH3

5 3 5 6 # CT1 HB

6 4 5 7 # CC CT1

7 5 5 8 # CT1 CT2

8 6 7 23 # CC OC

9 6 7 22 # CC OC

10 7 8 9 # CT2 HA

11 7 8 10 # CT2 HA

12 8 8 11 # CA CT2

13 9 11 12 # CA CA

14 9 11 20 # CA CA

15 10 12 13 # CA HP

16 9 12 14 # CA CA

17 10 14 15 # CA HP

18 9 14 16 # CA CA

19 10 16 17 # CA HP

20 9 16 18 # CA CA

21 10 18 19 # CA HP

22 9 18 20 # CA CA

23 10 20 21 # CA HP

head for the ch2lmp data.file:

LAMMPS (9 Jan 2009)

units real

atom_style full

boundary p p f

pair_style lj/charmm/coul/long 15.0 16.0 16.0

bond_style harmonic

angle_style charmm

dihedral_style charmm

improper_style harmonic ## this block is missing in the above version

kspace_style pppm 0.0001

kspace_modify slab 3.0

read_data phe_car.data # here is the name of the file above just sanning data file

3 = max bonds/atom # wrong output !!!

6 = max angles/atom

12 = max dihedrals/atom # wrong output

2 = max impropers/atom

orthogonal box = (0 -10 0) to (32 32 32.86)

data.file:

Bonds

1 5 1 3 # CC OC

2 5 1 2 # CC OC

3 4 1 8 # CC CT1

4 8 4 8 # CT1 NH3

5 10 5 4 # HC NH3

6 10 6 4 # HC NH3

7 10 7 4 # HC NH3

8 7 8 9 # CT1 HB

9 6 10 8 # CT1 CT2

10 9 10 11 # CT2 HA

11 9 10 12 # CT2 HA

12 2 13 10 # CA CT2

13 3 14 15 # CA HP

14 1 14 13 # CA CA

15 1 16 14 # CA CA

16 3 16 17 # CA HP

17 1 18 22 # CA CA

18 1 18 16 # CA CA

19 3 18 19 # CA HP

20 1 20 13 # CA CA

21 3 20 21 # CA HP

22 1 22 20 # CA CA

23 3 22 23 # CA HP

sjplimp · August 31, 2009, 1:49pm

But the discrepancy concerning the max bonds/atoms remains.

This value is important and is computed by LAMMPS by scanning your
data file and counting bonds on each atom. So if it is different
for 2 different runs using 2 different data files, it is certainly the
case that something in those data files is causing it. You mentioned
something about atom IDs being changed. If that's the case, then
it may be that the bonds assigned to a particular atom are different.

Steve

sjplimp · August 31, 2009, 1:50pm

One additional thought. If you have newton bond on (the default),
then a bond is assigned to one of the 2 atoms in it. If the
2 data files reversed the order of bonds then you could see a different
tally of which bonds are assigned to which atoms and thus the
count of max bonds per atom. But it wouldn't make a difference
since which atom a bond is assigned to is actually arbitrary.

Steve

_Sabine_Leroch · August 31, 2009, 2:36pm

Hi Steve.

Thank you for responding. I really carefully checked whether the bonds in both files are the same (the atomID are, as you say just names and can be assigned arbitrarily as long as both files give the same bonds per atom). Below there is an excerpt of the Bonds part of the data.file. I thought lammps (read.data.cpp) would scan the third and forth column and count how often a certain atomID is present. In this special case for example atomId 4 should be counted 4 times. So maxbonds/atom should be 4, is not it?

1 5 1 3 # CC OC

2 5 1 2 # CC OC

3 4 1 8 # CC CT1

4 8 4 8 # CT1 NH3

5 10 5 4 # HC NH3

6 10 6 4 # HC NH3

7 10 7 4 # HC NH3

8 7 8 9 # CT1 HB

9 6 10 8 # CT1 CT2

10 9 10 11 # CT2 HA

11 9 10 12 # CT2 HA

Best regards

Sabine

e*ve [email protected]> *8/31/2009 3:50 PM >>>
One additional thought. If you have newton bond on (the default),
then a bond is assigned to one of the 2 atoms in it. If the
2 data files reversed the order of bonds then you could see a different
tally of which bonds are assigned to which atoms and thus the
count of max bonds per atom. But it wouldn’t make a difference
since which atom a bond is assigned to is actually arbitrary.

Steve

_Sabine_Leroch · August 31, 2009, 3:50pm

By the way if I partly exchange columns in the list below i.e.

1 5 1 3 # CC OC

2 5 1 2 # CC OC

3 4 1 8 # CC CT1

4 8 4 8 # CT1 NH3

5 10 4 5 # HC NH3

6 10 4 6 # HC NH3

7 10 4 7 # HC NH3

8 7 8 9 # CT1 HB

9 6 10 8 # CT1 CT2

10 9 10 11 # CT2 HA

11 9 10 12 # CT2 HA

I receive max bonds/atom=4. Normally there should by no difference for both lists. There must be a bug.

Best regards

Sabine

“Sabine Leroch” <sabine.leroch@…1181…> 8/31/2009 4:36 PM >>>

Hi Steve.

Thank you for responding. I really carefully checked whether the bonds in both files are the same (the atomID are, as you say just names and can be assigned arbitrarily as long as both files give the same bonds per atom). Below there is an excerpt of the Bonds part of the data.file. I thought lammps (read.data.cpp) would scan the third and forth column and count how often a certain atomID is present. In this special case for example atomId 4 should be counted 4 times. So maxbonds/atom should be 4, is not it? But lammps counts only 3 bonds/atom!!!

1 5 1 3 # CC OC

2 5 1 2 # CC OC

3 4 1 8 # CC CT1

4 8 4 8 # CT1 NH3

5 10 5 4 # HC NH3

6 10 6 4 # HC NH3

7 10 7 4 # HC NH3

8 7 8 9 # CT1 HB

9 6 10 8 # CT1 CT2

10 9 10 11 # CT2 HA

11 9 10 12 # CT2 HA

Best regards

Sabine

e*ve [email protected]> *8/31/2009 3:50 PM >>>
One additional thought. If you have newton bond on (the default),
then a bond is assigned to one of the 2 atoms in it. If the
2 data files reversed the order of bonds then you could see a different
tally of which bonds are assigned to which atoms and thus the
count of max bonds per atom. But it wouldn’t make a difference
since which atom a bond is assigned to is actually arbitrary.

Steve

sjplimp · September 1, 2009, 2:23pm

Again, if you have newton bond on, which is the default, then
a bond is stored with only one atom, not with both atoms in the
bond. And it will be stored only with the first atom of the 2 listed.
So your examples all make sense to me; they are getting the
right answer.

Steve

_Sabine_Leroch · September 1, 2009, 3:06pm

Dear Steve,

ok thank you, I think I got it now, how lammps calculates maxbonds/atom. I saw it also by a printout I made today from the source code. One last thing I need to be sure of: I hope I understood your email from yesterday correctely. There you mentioned it has no influence on the results to which atom the bond is assigned. This means list 1 and 2 below lead to the same result despite the differences in maxbonds/atoms. Is this what you meant?

list 1:

1 8 4 8 # NH3 CT1> >

2 10 5 4 # HC NH3 > >

3 10 6 4 # HC NH3 > >

4 10 7 4 # HC NH3

list2:

1 8 4 8 # NH3 CT1> >

2 10 4 5 # NH3 HC> >

3 10 4 6 # NH3 HC> >

4 10 4 7 # NH3 HC

Best regards

Sabine

Steve Plimpton [email protected] 9/1/2009 4:23 PM >>>
Again, if you have newton bond on, which is the default, then
a bond is stored with only one atom, not with both atoms in the
bond. And it will be stored only with the first atom of the 2 listed.
So your examples all make sense to me; they are getting the
right answer.

Steve

sjplimp · September 2, 2009, 2:18pm

yes, which atom the bond is assigned to, or which
order the 2 atoms appear in the data file, makes
no differenence. There could be small round-off
distances that would accumulate in a long run,
but that's not an error

Steve