Problem of simulating gas adsorption in nanotube with lammps

I am simulating the adsorption of CH₄ and N₂ in a BN nanotube. CH₄ is modeled using a single-point LJ model, while N₂ follows the TraPPE force field, which consists of a 2LJ + 3-charge model. Specifically, CH₄ is represented as a single-point LJ model with one atom type, whereas N₂ is represented as a 2LJ + 3-charge model with two atom types.

I have a swnt.data file that contains the structure of the BN nanotube. This file includes only two atom types corresponding to B and N. A portion of its contents is as follows:

# LAMMPS data file converted from XYZ

        936 atoms

5 atom types

0.0        20.73053 xlo xhi
0.0        20.73053 ylo yhi
0.0        50.00000 zlo zhi

Masses

1 10.81
2 14.01
3 16.04  # CH4 (single-point LJ model)
4 14.01  # N atom in N2
5 1.0e-10  # Virtual site in N2 (mass = 0)

Atoms

    1   1   1   0.4000       16.58441        8.29221        1.20490
    2   2   1  -0.4000       16.45831        9.73271        1.20490
    3   1   1   0.4000       16.30191       10.43841        2.45730
    4   2   1  -0.4000       15.80731       11.79721        2.45730
    5   1   1   0.4000       15.47351       12.43831        1.20490
    6   2   1  -0.4000       14.64401       13.62281        1.20490
    7   1   1   0.4000       14.15571       14.15571        2.45730

The LAMMPS input script is as follows:

# --------------------- Initialization ---------------------
units real
atom_style full  # Use full to support charge and molecule ID
dimension 3
boundary s s p  # Non-periodic in x, y; periodic in z

# --------------------- Read BN Nanotube ---------------------
read_data swnt.data  # Directly load BN nanotube structure

# --------------------- Define Regions ---------------------
region nanotube cylinder z 10.365 10.365 8.29221 0 50 units box
region CH4_ads_region cylinder z 10.365 10.365 8.29221 0 50 units box
region N2_ads_region cylinder z 10.365 10.365 8.29221 0 50 units box

# --------------------- Force Field Parameters ---------------------
pair_style lj/cut/coul/cut 12.0  # 12 Å cutoff
pair_modify mix arithmetic

# No interaction between BN atoms
# Interactions between adsorbates and BN
pair_coeff 1 3 83.96 3.591  # B-CH4
pair_coeff 2 3 104.18 3.548 # N-CH4
pair_coeff 1 4 41.48 3.382  # B-N2
pair_coeff 2 4 51.21 3.338  # N-N2
pair_coeff 1 5 0.0 0.0      # B-Virtual Site
pair_coeff 2 5 0.0 0.0      # N-Virtual Site

# Adsorbate-adsorbate interactions
pair_coeff 3 3 148.0 3.73   # CH4-CH4
pair_coeff 4 4 36.0 3.31    # N-N (TraPPE)
pair_coeff 5 5 0.0 0.0      # Virtual site no interactions

# --------------------- BN Nanotube Charges ---------------------
set type 1 charge 0.4  # B atoms
set type 2 charge -0.4 # N atoms

# --------------------- Define Molecular Templates ---------------------
molecule CH4_mol CH4_molecule.txt
molecule N2_mol N2_molecule.txt

# --------------------- GCMC Adsorption Simulation ---------------------
fix gcmc_CH4 all gcmc 1000 10000 10000 3 29494 298.0 0.1 region CH4_ads_region fugacity 10.0 mol CH4_mol
fix gcmc_N2 all gcmc 1000 10000 10000 4 54321 298.0 0.1 region N2_ads_region fugacity 10.0 mol N2_mol

# --------------------- Simulation Control ---------------------
thermo 1000
thermo_style custom step temp press pe ke etotal density vol
dump mydump all atom 1000 dump.lammpstrj
run 500000

The CH4_molecule.txt file is as follows:

# CH4 molecule (single-point LJ model)

1 atoms

Atoms

1 3 0.0 0.0 0.0 0.0

The N2_molecule.txt file is as follows:

# TraPPE N2 molecule definition (corrected version)
3 atoms  # Now includes 2 N atoms + 1 virtual site
1 bond
3 charges

Atoms
1 4 -0.482  0.0  0.0 -0.55  # First N atom with partial charge
2 4 -0.482  0.0  0.0  0.55  # Second N atom with partial charge
3 5  0.964  0.0  0.0  0.0   # Virtual site with +0.964 e charge

Bonds
1 1 2  # N-N bond

Error Message

Upon running LAMMPS, I encountered the following error:

LAMMPS (29 Aug 2024 - Update 1)
Reading data file ...
  orthogonal box = (0 0 0) to (20.73053 20.73053 50)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  936 atoms
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     0 = max # of 1-2 neighbors
     0 = max # of 1-3 neighbors
     0 = max # of 1-4 neighbors
     1 = max # of special neighbors
  special bonds CPU = 0.000 seconds
  read_data CPU = 0.004 seconds
Setting atom values ...
  936 settings made for charge
Setting atom values ...
  0 settings made for charge
ERROR on proc 0: Unknown section 'Atoms' in molecule file
 (src/molecule.cpp:664)
Last command: molecule CH4_mol CH4_molecule.txt

First off, you posted to the wrong forum. This is a LAMMPS question, it should go into one of the LAMMPS categories. Since this is your first post here, I’ve corrected it during the approval process for new user posts. Don’t expect this to happen in the future. Posting to the wrong forum will likely have your post ignored.

Molecule files do not have an Atoms section, only data files. You cannot pass a data file as a molecule file. While they look similar, their format is different.

Thanks. Please see that the molecule files have Atoms section.

You misunderstood what I am saying. I wasn’t saying that they should have one, but they should not. Please see the documentation for the molecule command.

As I said, you seem to be trying to pass off a data file as a molecule file and that will not work and will get you exactly the error that you are seeing.

Another problem will be that both, molecule and data files must have an empty line following the header. Thus your files will have that line ignored and thus LAMMPS will complain about missing lines.

Why should not have a Atom section in CH4_molecule.txt?

Because that is the documented behavior and that is why you get an error.

Can you modify the CH4_molecule.txt and N2_molecule.txt? thanks.

Please read the documentation for the molecule command and make the changes yourself.

It is not my job to do your work.

I read the doc and solve the CH4_molecule.txt error. But I got the error for N2_molecule.txt.

TraPPE N2 molecule definition (corrected version)

3 atoms # Now includes 2 N atoms + 1 virtual site
1 bonds

Coords

1 0.0 0.0 -0.55 # First N atom with partial charge
2 0.0 0.0 0.55 # Second N atom with partial charge
3 0.0 0.0 0.0 # Virtual site with +0.964 e charge

Types

1 4 # First N atom
2 4 # Second N atom
3 5 # Virtual site

Charges

1 -0.482
2 -0.482
3 0.964

Bonds

1 1 1 2

The error is below:
LAMMPS (29 Aug 2024 - Update 1)
Reading data file …
orthogonal box = (0 0 0) to (20.73053 20.73053 50)
1 by 1 by 1 MPI processor grid
reading atoms …
936 atoms
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.011 seconds
Setting atom values …
936 settings made for charge
Setting atom values …
0 settings made for charge
Read molecule template CH4_mol:

CH4 molecule (single-point LJ model)

1 molecules
0 fragments
1 atoms with max type 3
0 bonds with max type 0
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
ERROR: Invalid bond type in Bonds section of molecule file: 1 1 1 2 (src/molecule.cpp:1087)
Last command: molecule N2_mol N2_molecule.txt

I can not understand the second number in the line

Bonds

1 1 1 2

Isn’t it 1??

Yes, but you cannot define a bond when you have not reserved any space for storing bonds. That has to be done when the simulation box is defined and in your case that is with the read_data command. There are the various “extra” parameters for that command. Again, this is documented behavior.