Inserting molecules into a LAMMPS (output) data file with constraints?

LAMMPS
1

Dear LAMMPS Users and Developers,

I have a LAMMPS input file that creates two copper slabs (below each other) and then minimize the structure by cg and after that I run an NVT simulation using Nosé-Hoover thermostat.

My goal is to insert oxygen molecules into this optimized structure in a specific way:
-oxygen molecules should be at least 3.5 Å from each other
-oxygen molecules should be at least 6 Å from copper atoms

I checked the molecule command but do not know how I could define the distances between the molecules-molecules (O2-O2) and molecules-atoms (Cu-O2).

I use ReaxFF.

Thank you for your help,
Patrik

PS. I could not attach my file so I copied the input:

LAMMPS Input Script for Two Flat Copper Slabs with ReaxFF Potential and Charge Equilibration

Set units and atom style

units real
atom_style charge

Define the simulation box

lattice fcc 3.615 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region total_box block -161 161 -161 161 -40.5 40.5 units lattice
create_box 1 total_box

Create the first copper slab in the lower half of the box

region slab1 block -20 -8 -20 -8 -1.5 1.5 units lattice
create_atoms 1 region slab1
mass 1 63.546

Create the second copper slab in the upper half of the box

region slab2 block -20 -8 -20 -8 -6 -3 units lattice
create_atoms 1 region slab2

Define the reactive force field

pair_style reaxff lmp_control
pair_coeff * * ffield.reax.cuo Cu

Charge equilibration for ReaxFF

fix 2 all qeq/reax 1 0.0 10.0 1e-6 param.qeq

Minimize the structure using conjugate gradient method

min_style cg
minimize 1e-4 1e-6 1000 10000

Define the simulation settings

timestep 0.001
thermo 100

Apply Nosé-Hoover thermostat for NVT ensemble to equilibrate at 300K

fix 1 all nvt temp 300.0 300.0 50

Run the equilibration

run 10000

Output final configuration

write_data final_configuration_middle.data

2

You can try using the create_atoms command with:

  • a molecule template (to insert molecules)
  • the overlap keyword with a value of 3.5 A (to avoid overlap between atoms)
  • the region style (to place your molecule only in a region far from the walls)

edit: sorry I meant the random style with a region-ID, not the region style

1 Like
3

On top of @simongravelle suggestions:

  • Define a region that is exactly 3.5 A above the Cu surface.
  • From the volume of the region to be filled with O2 molecules and the density you want to achieve, determine how many molecules to insert.
  • You can’t attach files, but you can read the forum guidelines. Please format your post accordingly.
1 Like
4

Sorry, how should this command look like?
Can the oxygen be created by molecule command if I would like to apply charge atom_style?

I mean I get an error (" Invalid bond type in Bonds section of molecule file: 1 1 1 2 (src/molecule.cpp:1087) ") for the Bonds style for the following oxygen molecule:

#LAMMPS oxygen DATA file
2 atoms
1 bonds
0 angles

Coords

1 -1.49100 2.13600 0.01500  # ID Type Charge x y z
2 -0.29300 1.59000 0.01500  # ID Type Charge x y z

Types

1	1 #O
2	1 #O

Charges

1	-0.82000
2	-0.82000


Bonds

1 1 1 2

Angles
5

For ReaxFF you do not define bonds since ReaxFF computes them on the fly. So your molecule file also must not define a bond yet your header says there is one bond.

1 Like
6

Thank you for your answer!

If I do not define Bonds then I have to load the parameters for oxygen in the line pair_coeff. When I try it, I get the error:

ERROR: Incorrect args for pair coefficients (src/REAXFF/pair_reaxff.cpp:280)
Last command: pair_coeff * * ffield.reax.cuo Cu O # placeholder, change with actual coefficients

I used the following settings:

pair_style reaxff lmp_control  # change as per your system
pair_coeff      * * ffield.reax.cuo Cu O # placeholder, change with actual coefficients

The force field file contains Cu and O exactly with these characters.
What did I do wrong?

7

It’s difficult to analyze errors from small pieces of code, but if I refer to your original post, it seems that you have only 1 atom type, but the way you write your pair_coeff suggests you have 2. Therefore, LAMMPS is complaining.

Note that this issue is unrelated to your original question, and therefore should not be part of the same post.