How to add atom attribute molecule for each layer in LAMMPS?

Question

I want to model several layers of hBN in LAMMPS.

First, for two layers, I created a lammps data file from vasp POSCAR using atomsk. It is like this:

 # hBN_bilayer
  
        3200  atoms
           2  atom types
 
      0.000000000000      50.248565674000  xlo xhi
      0.000000000000      87.033068756000  ylo yhi
      0.000000000000      24.853632450000  zlo zhi
 
Masses
 
            1   10.81000000             # B
            2   14.00700000             # N
 
Atoms # atomic
 
         1    2        0.000000000000       1.450551189000       5.780448675000
         2    2        1.256214141900       3.626377907900       5.780448675000
         3    1        0.000000000000       1.450551189000       1.926816225000
...

I want to use extep potential inside each layer and ilp/graphene/hbn between layers.

pair_style hybrid/overlay extep ilp/graphene/hbn 16.0 1
pair_coeff * * extep BN.extep B N
pair_coeff * * ilp/graphene/hbn BNCH.ILP B N

But when I run this, it says

ERROR: Pair style ilp/graphene/hbn requires atom attribute molecule

How do I specify this attribute for different layers?

I tried to read the manual page https://docs.lammps.org/pair_ilp_graphene_hbn.html but I could not get how to do it.

Answer 1

The 'molecule' identifier will need be a separate column field for each atom. In your case, each layer will count for one molecule. So layer 1 will have a 1 in that column and so on.

From the documentation of your pair_style: "To keep the intralayer properties unaffected, the interlayer interaction within the same layers should be avoided. Hence, each atom has to have a layer identifier such that atoms residing on the same layer interact via the appropriate intralayer potential and atoms residing on different layers interact via the ILP. Here, the molecule id is chosen as the layer identifier, thus a data file with the 'full' atom style is required to use this potential."

One way of writing that column into the lammps input file is with the python package ASE's ase.io.write() method, see for example https://wiki.fysik.dtu.dk/ase/ase/io/formatoptions.html#lammps-data

In particular, you will want to specify atom_style='full' when writing with ASE. You should get something like this with the "full" format:

Atoms # full
 
   1           1          2        0.00 0.000000       1.450551       5.780448
   2           1          2        0.00 1.256214       3.626377       5.780448
   3           2          1        0.00 0.000000       1.450551       1.926816
# id molecule-id  atom-type           q        x              y              z

The extra columns are the index of molecules (molecule-id) or layers to which each atom belongs, plus the charge q. You should check that it makes sense before starting the simulation.

Answer 2

I accepted Andrey Poletayev's answer and upvoted it because it gives me the right direction, although it needs to be corrected (see my comment).

However, I couldn't do it with the python package ASE, so I wrote a function in Julia that converts an atomic style LAMMPS data file to a full style.

Because my question was upvoted, I think other people might be interested in it, so I post my code here

function atomic2full(input_atomic, output_full)
    lines = open(input_atomic) do file; readlines(file); end

    # Get the number of the line "Atoms # atomic"
    idx = findfirst(x -> startswith(x, "Atoms"), lines)

    # Number of atoms
    n_atoms = parse(Int, filter(x -> '0' <= x <= '9', lines[3]))

    # Indices of the first and the last lines with atomic coordinates
    ib = idx + 2; ie = idx + 2 + (n_atoms - 1);

    # Array of z coordinates
    z_coords = Array{Float64}(undef, n_atoms)
    z_coord_position = 5
    for (iz, iline) ∈ enumerate(ib:ie)
        z_coords[iz] = parse(Float64, (split(lines[iline])[z_coord_position]))
    end

    # Unique z values up to two digits after the decimal point
    z_uniques = unique(trunc.(z_coords, digits = 2))
    sort!(z_uniques)

    q = "0.0" # charge

    for (iz, iline) ∈ enumerate(ib:ie)
        molecule_tag = findfirst(x -> isapprox(x, z_coords[iz]; atol=0.01), z_uniques) 
        line_vec = split(lines[iline])
        insert!(line_vec, 3, q)
        insert!(line_vec, 2, string(molecule_tag))
        lines[iline] = "      "*join(line_vec, "      ")
    end

    lines[idx] = replace(lines[idx], "atomic" => "full")
    
    open(output_full, "w") do file
        for line ∈ lines
            write(file, line*"\n") 
        end
    end    
end;  

Usage:

atomic2full(input_filename_atomic, output_filename_full)