Classical rigid LJ: Non-numeric atom coords

LAMMPS LAMMPS General Discussion
,
1

Hi everyone:

I want to politely ask a question about applying LAMMPS simulation to a rigid LJ molecule with an NPT ensemble:

Here’s the input file:

INPUT

units real
atom_style full
boundary p p p
pair_style lj/cut/coul/cut 15.0
bond_style harmonic
pair_modify mix geometric tail yes

read_data DATA.data
bond_coeff * 1050.000000 1.102000
group clump union all
timestep 0.5
fix 1 clump rigid/npt molecule temp 64.000000 64.000000 5.0 iso 0.493462 0.493462 500.0

dump 1 all atom 100000 INPUT.1.dump #0.1ns per dump
dump_modify 1 scale no image yes
thermo_style custom step time cpu tpcpu etotal ke pe ebond eangle edihed eimp evdwl ecoul elong enthalpy temp press vol density
thermo 5000
thermo_modify flush yes
restart 1000 INPUT.1.rst INPUT.2.rst

run 10000000 # 5ns
write_restart INPUT.rst

With the data file that results in the error:

LAMMPS data, DATA file1 with filename DATA.data

2 atoms
1 bonds
0 angles
0 dihedrals
0 impropers

1 atom types
1 bond types
0 angle types
0 dihedral types
0 improper types

0.000000 500.000000 xlo xhi
0.000000 500.000000 ylo yhi
0.000000 500.000000 zlo zhi

Masses

1 14.006700

Atoms

1 1 1 0 193.593092 158.841510 31.234906
2 1 1 0 194.593659 159.193643 30.696596

Bonds

1 1 2 1

Pair Coeffs

1 0.075900 3.262100

Bond Coeffs

1 1050 0.969900

This should be a straightforward system, and the two atoms are not overlapping, yet lammps is resulting in the following error:

ERROR on proc 0: Non-numeric atom coords - simulation unstable (src/OPENMP/domain_omp.cpp:58)
Last command: run 10000000 # 5ns

The input files should be correct because when shifting this single molecule in the x-axis, the issue magically disappears for some reason. What could went wrong with the previous data file?

The following data file will somehow not result in an error:

LAMMPS data, DATA file2 with filename DATA.data

2 atoms
1 bonds
0 angles
0 dihedrals
0 impropers

1 atom types
1 bond types
0 angle types
0 dihedral types
0 improper types

0.000000 500.000000 xlo xhi
0.000000 500.000000 ylo yhi
0.000000 500.000000 zlo zhi

Masses

1 14.006700

Atoms

1 1 1 0 93.593092 158.841510 31.234906
2 1 1 0 94.593659 159.193643 30.696596

Bonds

1 1 2 1

Pair Coeffs

1 0.075900 3.262100

Bond Coeffs

1 1050 0.969900

When applying the same forcefield in a system with more molecules, this error persists. Could there be a format error in the data file? Any help or suggestion is highly appreciated!!

2

Here’s an example run if anyone wants to have a try:

https://drive.google.com/drive/folders/1oE6T_TaAlV8RmpYqRZnP8695yTgpxRG0?usp=sharing

3

What if you use either fix rigid/nvt or fix rigid/nve? I am very suspicious of the use of a barostat (or thermostat!) with a very large box containing a single rigid object.

4

Hi @oliver_song,

Some points that might help clarify this issue:

  • The error message refers to the file domain_omp.cpp and the line involved. You can take a look at the code (or the similar function in domain.cpp to get what’s going on and what check the code does before throwing the error. Also take a look at the assumptions on barostating algorithms and how pressure is defined in periodic systems. You will realise why:
    1. It does not make much sense to compute pressure and equilibrate a box with a single molecule inside it. What will be the pressure of a single molecule that interacts with nothing? How is pressure expected to be relaxed? In what aspect does moving a molecule might change anything?
    2. Anything can numerically happen in these cases and, well, sometime it will work, sometime it won’t, but as the algorithms are not made for these kind of systems, this is not very telling on what to expect on real productions simulations.
  • If single molecule systems are trivial systems (depends on what you expect but let’s assume here that it is the case), barostat and thermostat algorithms are anything but trivial. And they’re built up from less to more sophisticated NVE → NVT → NPT. This is a crude description, but if you try to nail a bug or an issue, “simplify simplify simplify”, start simple and try to understand what’s going on at each step and what is the difference with the next one that might make anything go wrong. See @srtee’s comment (BTW, no integration issues in both other cases but I only tested one of the two single molecule systems).
  • If I understand things correctly you went from the case you provide in the google drive to single molecule systems. I think you’re on the wrong path. The reason I think your system is throwing an error in rigid/npt is because you start with a low density gas phase system with zero kinetic energy which is very far away from your expected phase space (a guess from the parameter of your system but the system goes to higher density anyway). Basically this is a bad starting point (90% of the reason LAMMPS throws this error from my own forum post reading). As mentioned this issue has been discussed extensively on the forum and the mailing list. However it is more about system preparation and physics. and you will get more significant insight through direct discussion and comments on what you try to achieve with your advisor or colleagues.
5

Hello @Germain and @srtee

I want to show my tremendous appreciation for the suggestions from both of you.

Yeah, you’re right. The use of the barostat or thermostat likely caused the problem. When I reduce the timestep to 0.1 (in my understanding), it seems like the coordinates of the atoms are not aggressively changed by them, and the error is no longer there.